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Abstract 

We evaluate numerically and analytically the dynamic critical exponent z for five gauge- 
fixing algorithms in SU{2) lattice Landau-gauge theory by considering the case (3 = oo. 
Q\ , Numerical data are obtained in two, three and four dimensions. Results are in agreement 

with those obtained previously at finite (3 in two dimensions. The theoretical analysis, valid 
for any dimension d, helps us clarify the tuning of these algorithms. We also study gen- 
eralizations of the overrelaxation algorithm and of the stochastic overrelaxation algorithm 
and verify that we cannot have a dynamic critical exponent z smaller than 1 with these 
local algorithms. Finally, the analytic approach is applied to the so-called A- gauges, again 
at (3 = oo, and verified numerically for the two-dimensional case. 
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1 Introduction 



Lattice gauge fixing is a necessary step in our understanding of the relationship between con- 
tinuum and lattice models. In fact, the continuum limit is the weak-coupling limit and a weak- 
coupling expansion requires gauge fixing. Thus, even though gauge fixing is in principle not 
needed for the lattice formulation of QCD, one is led to consider gauge-dependent quantities on 
the lattice as well, such as gluon, ghost and quark propagators, vertices, etc. [1]. It is therefore 
important to devise numerical algorithms to gauge-fix a lattice configuration efficiently. The 
main issue regarding the efficiency of these algorithms is the problem of critical slowing- down 
(CSD), which occurs when the relaxation time of an algorithm diverges as the lattice volume 
is increased (see for example [2, 3]). Besides the problem of CSD, we are also interested in un- 
derstanding which quantities should be used to test the convergence of the gauge fixing and in 
finding prescriptions for the tuning of parameters in the algorithms, when tuning is needed. 

In Ref. [4] we have studied the problem of critical slowing-down for five gauge-fixing algo- 
rithms in SU (2) lattice Landau-gauge theory on two-dimensional lattices with periodic boundary 
conditions. We obtained that the local method called Los Alamos has dynamic critical exponent 
z pa 2, the three improved local methods we considered — the overrelaxation method, the stochas- 
tic overrelaxation method and the so-called Cornell method — have critical exponent z«l, and 
the global method of Fourier acceleration completely eliminates critical slowing-down. All these 
methods, except for the Los Alamos method, involve tuning. In that reference we also reported 
a detailed discussion and analysis of the tuning for the overrelaxation, the stochastic overrelax- 
ation and the Cornell methods, and we made a comparative study of several quantities usually 
employed in the literature to check the convergence of the gauge fixing. 

Here we redo that analysis for the case /3 — oo, which can be studied analytically, and we 
include test runs in two, three and four dimensions 1 . Let us recall the numerical problem we want 
to study: for a given (i.e. fixed) lattice configuration {U^(x)}, we look for a gauge transformation 
{g(x)} that brings the functional 

£u[g] - (i-i) 

U V fj,= l X z 

= + (1.2) 

= l-^EEll^^ + ^ffW] (1-3) 

a V ft=l X Z 

to a local minimum, starting from a randomly chosen configuration {g(x)}. Here, x (with coor- 
dinates — 1,2, . . . , N) are sites on a <i-dimensional lattice with periodic boundary conditions, 
V = N d is the lattice volume 2 and U^(x) and g(x) are SU(2) matrices. In order to analyze 
the CSD of an algorithm we have to investigate if, and with what exponent, its relaxation time 
t diverges as the lattice size increases. To this end, we have to evaluate r for different lattice 
sides N at "constant physics", namely as the lattice size is increased, the physical size of the 
lattice should remain fixed. This is done by introducing a correlation length £ and by keeping 
the ratio iV/£ constant. The lattice configuration {U^x)} is usually chosen in equilibrium with 

1 Partial results can be found in [5] . 

2 In order to simplify the notation we don't consider asymmetric lattices. 
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the Gibbs weight exp (-S \J3, £/]), where S \J3, U] is the standard Wilson action in d dimensions. 
Since S [(3, U] depends only on the coupling (3 and on the lattice side N, we must change (3 
and N in such a way that the pairs (/3,N) correspond to the same value of N/£. For example 
- in two dimensions and for £ = l/y/it, where k is the string tension — £(/3, TV) is known in 
the infinite- volume limit [6] and for large (3 we have the simple relation oo) oc (3 1 ^ 2 . [In 
Ref. [4] we considered different pairs ((3,N) with the ratio N 2 / f3 fixed and equal to 32, which 
corresponds to iV/£ rs 7.] Of course, if the function N) is not available, one should evaluate 
directly (i.e. numerically) the correlation length £, and tune (3 and N so that the ratio iV/£ is 
kept (approximately) fixed. 

Since we are interested in studying gauge-fixing algorithms, i.e. minimizing the functional 
(1.2) for a given lattice configuration {U^(x)}, a simpler possibility [7] is to work with lattices 
at P — oo. This corresponds to using the Gibbs weight U Xjtl 5(1. — U^(x)), namely the lattice 
configuration is completely ordered. In this case, the string tension k is zero for any lattice side 
N, and therefore N/£ = N^/k is constant and equal to zero for all the pairs ((3 = oo, N). Note 
that with this particular choice the link variables U^x) are set to 1L and the configuration {g(x)} 
is chosen randomly. Thus, Ujf\x) is a pure gauge configuration, i.e. U^\x) = g(x) g~ x (x + e M ), 
which should be driven by the gauge fixing to the the vacuum configuration Ujf\x) = 1L. As we 
will see in Section 3, the advantage of using this particular case — for which we know the final 
gauge-fixed configuration — is that we can study analytically the efficiency of the gauge-fixing 
algorithms. 

In Section 2 we briefly review the five gauge-fixing algorithms considered in Ref. [4]. The 
case (3 = oo is studied analytically for general dimension d in Sections 3, 4 and 5 and the 
results are checked numerically in two, three and in four dimensions in Section 6. The analytic 
approach (see Section 3) is done by mapping — in the limit of large number of gauge-fixing 
sweeps t — the original problem (1.2) into the minimization of the action of a three- vector 
massless-free-field model [see eq. (3.6)]. In this way we obtain that the Cornell method coincides 
with the overrelaxation method, in agreement with the result found in Ref. [4], and we show 
that the Fourier acceleration method can eliminate critical slowing-down completely. In Section 
4 we review the analysis of CSD done in Refs. [8, 9] for the thermalization of the Gaussian 
model and we apply it to the four local gauge-fixing algorithms considered here. This analysis 
confirms the results obtained in Ref. [4] and will give us a better understanding of the tuning 
for the three improved local algorithms. Moreover, in Section 5, we study generalizations of 
the overrelaxation algorithm and of the stochastic overrelaxation algorithm and check that we 
cannot have a dynamic critical exponent z smaller than 1 with these local algorithms. Note that 
this result applies directly to the problem of thermalizing a massless free field [8, 9, 10]. Finally, 
the analytic approach is extended to the so-called A gauges (see for example [11]) in Section 7 
and the results are verified numerically in the two-dimensional case. 

Let us remark that the case /3 = oo gives information valid also for the algorithms at finite 
f3. In fact, for each algorithm and "constant physics", the evaluated relaxation times should be 
fitted by a function of the type 

r = cN z , (1.4) 

where the dynamic critical exponent z should not depend on the constant physics, i.e. it should 
be the same 3 at finite f3 and at f3 — oo. On the contrary, the constant c should depend on the 

3 Of course, if the system undergoes a phase transition going from /3 = to [3 = oo, then an algorithm can 
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ratio N/£. In particular, for each algorithm, we should find a faster convergence of the gauge- 
fixing procedure — and therefore a smaller constant c — when the configuration U^(x) is less 
disordered, i.e. for larger values of (3. This implies 



c{(3 = oo) < c( finite (3) 



(1.5) 



In Section 8 we check this inequality in two dimensions using the data obtained in Ref. [4] at 
finite (3 and the data produced for the present work. In the same section we also give our final 
comments and conclusions. 



2 Gauge-Fixing Algorithms 

In this Section we review briefly the five gauge-fixing algorithms we want to study, i.e. the so- 
called Los Alamos method, the overrelaxation method, the stochastic overrelaxation method, the 
so-called Cornell method and the Fourier acceleration method. To this end, let us recall that 
the updated gauge transformation at each step of the algorithm can be written for the four local 
methods as follows 4 



g( LosAL \x) 



g { - cornell \x) 



EE h\x) 

aU{x)h\x) + [1 - ajV(x)T(x)/2] g(x) 



9 



(over) 



(stoc) 



(X) EE 



(X) 



^l + a 2 jV 2 (x)[l-T 2 (x)/4] 

uj w{x) + (1 — uj) g(x) 
y/l+u(u>-l)[2-T(x)] 

h^(x)T(x) — g(x) with probability p 
h){x) with probability 1 — p 



(2.1) 
(2.2) 

(2.3) 
(2.4) 



Notice that, with the exception of the Los Alamos method, all the above algorithms have a 
tuning parameter, i.e. a for the Cornell method, uj for the overrelaxation algorithm and p for the 
stochastic overrelaxation algorithm. Also, in eqs. (2.1)-(2.4) we have used the definition 



a 

K x ) = [ u »( x ) 9\ x + + u l( x - 9\ x - 



(2.5) 



and the fact that the matrix h(x) is proportional to an SU (2) matrix, namely it can be written 



as 



h(x) ee Af(x) h(x) 



with h(x) e SU(2) and 



jV(x) = yjdet h{x) = J—h(x)M(x) 



(2.6) 
(2.7) 



show a different behavior in different phases (see for example Ref. [12]). 
4 See Ref. [4] for more details. 
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We also define 
and 



T(x) 



w(x) = g(x) h(x) 
1 



X 



Trw(x) . 



Using the definition of Ujf\x) [see eq. (1.2)] and eq. (2.5) one can check that 



w(x) = £ U^(x)+uf(x-e,) 



(2.8) 
(2.9) 

(2.10) 



For the Fourier acceleration algorithm one usually writes [4] 
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(Fourier) . 



x) OC < 1L — a 



(2.11) 



where a is a tuning parameter, 1L is the identity matrix, F indicates the Fourier transform, F 1 
is its inverse, 



( V ■ A®) (x) = J2 Af{x) - A<*\x - e„) 



is the lattice divergence of the (gauge-transformed) gauge field 

1 

2 

and 



4 fl) (*) = \\u^\x)-uf{x) 



,,=i W / M= i 



(2.12) 



(2.13) 



(2.14) 



is the squared magnitude of the lattice momentum. The vector k has components (// = 
1, . . . , d) given by k^ N = 0, 1, . . . , N — 1 , where N is the lattice side. Sometimes eq. (2.11) is 
also written with p^ax / P 2 {k) instead of 1 / p 2 {k) . However, p^ a;E = 4rf is a constant and can 
therefore be absorbed into the tuning parameter a. 
Clearly, using eqs. (2.10), (2.12) and (2.13), we have 



(V-A^)(x) = lT,U^(x)-U^\x)-U^(x-e,) + uf(x-e,) (2.15) 

(2.16) 



/i=i 

w(x) — w^(x) 



For a matrix g G SU (2) we use the parametrization 

g = g 1L + ia ■ g = I 



go + m 92 + ig\ 
92 + igi go - m 



(2.17) 



where the components of a = (a 1 , a 2 , a 3 ) are the three Pauli matrices. Since w(x) is proportional 
to an SU(2) matrix we can write w(x) = w (x)lL + ia ■ w(x) and 



(V- A {9) ) (x) = ia-w(x) . 
5 



(2.18) 



Defining now 



we have 



u(x) = IF' 1 



p 2 (k) 



Fw 



g {Fourier \x) oc {1L - iaa-u(x)} g(x) 
and the Fourier acceleration update can be written as 

1L — i a a • u(x) 
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(Fourier) , 



x) 



1 + a 2 | u(x) 



■.g{x) . 



(2.19) 
(2.20) 

(2.21) 



In order to simplify the expression for u(x) in eq. (2.19), recall that (minus) the lattice Laplacian 
A [defined in eq. (3.7) below] is diagonal in momentum space, with elements given by p 2 (k). 
Therefore we have 

F' 1 -t^ttt F = (-A)' 1 (2.22) 



p 2 (k) 



and 



u(x) 



-A 



w 



(x). 



(2.23) 



Clearly, the inverse of the (minus) Laplacian is defined only in the sub-space orthogonal to the 
constant mode k = 0. From the previous equation it is evident that the Fourier acceleration 
method is actually a Laplacian preconditioning algorithm. 



3 The Case (3 = oo 

As said in the Introduction, we consider here the case in which all the link variables U^(x) are 
set equal to the identity matrix. Then, the minimizing functional (1.2) becomes 

1 d Tr r 
£ u[g] = 1 -J7tJ2J2^T g(x)g\x + e^) 

av H=l x z 

= 7/E y ^ E {[*(*) " 9(x + e,)} [g(x) - g(x + e^} (3.1) 
= 7/ E 73 E II 9(x) - g(x + e M ) f . (3.2) 

V X Z(1 fj,= l 

Here we have used Tr g = Trg^,gg^ = 1L and [using eq. (2.17)] 

^Tr (g g ] ) =g g + g.g = g-g = \\g\\ 2 , (3.3) 

where g = (go,g) is a four-dimensional unit vector. 

Equation (3.2) is very similar to the action of a four- vector massless free field g(x). The only 
difference is that our field must satisfy the constraint ||g(x)|| 2 = 1, namely we are dealing with 
an 0(4) nonlinear a-model. However, after a few sweeps, we expect the (gauge-transformed) 
gauge configuration Ujf\x) = g(x) g\x + e M ) to be very close to the vacuum configuration, 
i.e. we expect the matrices g(x) to be very close to a constant matrix. Since the minimizing 
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functional (1.3) is invariant under global gauge transformations, i.e. transformations in which 
g(x) is ^-independent, we can set this constant matrix equal to the identity matrix 1L. Therefore, 
in the limit of large number of gauge-fixing sweeps t we can write 5 



g(x) = 1+ iea-f(x) + 0{e 2 ) , 



(3.4) 



with e C 1. This approximation is justified because it is exactly in the limit of large t — in 
which only the slowest mode survives — that we evaluate the relaxation time of each algorithm. 
Now, by using eq. (3.4), we can rewrite the minimizing functional (3.2) as 



0(e 3 ) 



(3.5) 



and, using the periodicity of the lattice, 

*[/] = /(*)•(-*/)(*) +0(e 3 ), 

where —A is (minus) the lattice Laplacian, defined by 



(3.6) 



(-A/)(x) = £[2/(:r) - f(x + e„) - f(x - e„) 



(3.7) 



Therefore, by consistently keeping only terms up to order e in g(x), we have that £ [/] - at the 
lowest order in e — is the action of the three- vector massless free field f(x). In particular, we 
can update this field, in order to minimize the action, without taking into account the problem 
of the unitarity of g. In fact, from eq. (3.4) we have g(x) g\x) = 1L + 0(e 2 ) for any field f(x). 
It is not difficult to translate the update of the g(x) variables for the five algorithms considered 

— * 

here into an update for the f(x) field. To this end, we use the approximation in eq. (3.4) and 
obtain [see eqs. (2.5) and (2.8)] 



K x ) = E [dKx + e^) +g\x-e tl ) 

d 

= 2d! - iea- [fa + ej + fix-ej 
w(x) = g(x) h(x) 



+ 0{e 2 ) 



(3.8) 



2dl. + tea- J2[2f(x) - f(x + e M ) - f(x - e M ) 



+ 0{t 2 ) 



and 



w(x) 



e (-Af)(x) + 0(e 2 ) , 



(3.9) 



(3.10) 



5 Note that with the parametrization g(x) = go(x)l. + i e f(x) ■ a + 0(e 2 ) the unitarity condition for g{x) 
implies that the corrections to go{x) « 1 start at order e 2 . 
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where we used eq. (3.7). From these expressions it follows that [see eqs. (2.7), (2.9) and (2.23)] 



Af(x) = 2d + 0(e 2 ), 
T(x) = 2 + 0{e 2 ) , 
u(x) = ef(x) + C(e 2 ) . 
Actually, to be precise [see comment after eq. (2.23)], we should write eq. (3.13) as 

u(x) = e [f(x) -/o] + 0(e 2 ), 

where 



(3.11) 
(3.12) 
(3.13) 

(3.14) 
(3.15) 



However, the zero mode of g(x), and therefore of f(x), does not contribute to the value of the 
minimizing functional (3.1). Thus, in eq. (3.4), we can always consider a field f(x) with zero 
constant mode f . 

Now, by substituting eq. (3.4) and eqs. (3.8)-(3.13) into eqs. (2.1)-(2.4) and (2.21) we obtain 

— * 

for the updated f{x) variables (at the lowest order in e) 



7 (LosAl.) 



-*{cornell) 

f (x) 

7 (over) 



-i(stoc) 

f {x) 



-*(Fourier) 

/ {x) 



^J2[f(x + e,) + f(x-e,) 



2daf LosM \x) + (1 - 2 da) f(x) 



7 (LosAl.) 

uf (x) 



u) f(x) 



2 ^ (x) — f(x) with probability p 



7 (LosAl.) 

f (x) 



a) f(x) 



with probability 1 — p 



(3.16) 

(3.17) 
(3.18) 

(3.19) 
(3.20) 



The interpretation of these updates is clear if we consider the local minimization of £ [/], i.e. 
we consider all /(x)'s fixed for x ^ y, and we make explicit the dependence of the minimizing 
functional (3.5) on the value of the field / at site y, namely (to order e 2 ) 

^2 



£[/G/): 



V 



f(y) 



M - 2f 0sM \y) 



+ constant terms , 



(3.21) 



where we have used eq. (3.16). Since / 
square and write 

J2 



(LosAl.) 



(y) does not depend on f(y) we can complete the 



£[f(y)\ = y 



f(y) - f 



7 (LosAl.) 



(y) 



+ constant terms . 



(3.22) 



Then it is clear that the Los Alamos update brings this local functional to its minimum, while 

->(LosAl ) -* 

the choice f(y) — > 2 / (y) — f(y) does not change the value of £[/]. From formulae 
(3.16)-(3.20) it is also evident that: 
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• The Los Alamos method corresponds to the usual Gauss-Seidel method [13, 14]. 

• The Cornell method coincides with the overrelaxation method if we use the relation [see 
eqs. (3.11), (3.17) and (3.18)] 

oj = 2da = aAf(x) + C(e 2 ) . (3.23) 

This confirms the result found analytically and numerically in two dimensions at finite (5 
(see Sections 5 and 7.3 in Ref. [4]). 

• The overrelaxation method corresponds to the usual successive overrelaxation method [13, 
14]. 

• Since the vacuum configuration {U^{x) = 1L} corresponds to f(x) = for all x, it is clear 
that the Fourier method can minimize S [f] in only one step if we set a = 1. Thus, in this 
case, critical slowing-down is completely eliminated. Note that the tuning condition, i.e. 
a — 1, does not depend on the dimension d of the lattice. 

Moreover, it is now possible to study analytically the critical slowing-down for the four local 
algorithms following the analyses in Refs. [8, 9] (see next section). As we will see, the results of 
this analytic approach confirm the dynamic critical exponents found numerically in two dimen- 
sions at finite f3 [4]. Also, these results are particularly interesting with respect to the problem 
of tuning the improved local algorithms (see Section 6.2). 

Let us notice that for t going to infinity all components of f(x) must go to zero for the 
algorithm to converge. Therefore, if in some basis we can write 

/* = C/t-i , (3-24) 
where C is the updating matrix, then we should have 

lim || /J = lim HCVoll = 0, (3.25) 

t— >oo t-^oo 

for any reasonable definition of the norm ||/|| and of the initial condition Jq. From the definition 
of the norm of a matrix (see for example Section 1.5 in [14]) it follows that 

||C*/ || < W&W H/oll • (3.26) 

Thus, the limit in eq. (3.25) is verified if ||C*|| goes to zero when t goes to infinity, i.e. if the 
matrix C* goes to in the same limit. This happens if and only if (see theorem 1.4 in [14]) the 
spectral radius p(C) of the matrix C is smaller than 1, where 

p(C) = max |A| (3.27) 

Afc(T(C ) 

and o~(C) is the set of all the eigenvalues of the matrix C. One can also prove (see theorem 1.6 
in [14]) that 

p(C) = lim IICl 1 /* (3.28) 

for any matrix norm || . ||. It follows that the algorithm converges if all the eigenvalues of the 
updating matrix C are smaller than 1 (in absolute value), i.e. if p(C) < 1. From now on, we 
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consider only the case p{C) < 1. Then, we can define the relaxation time r > through the 
relation 

e- 1 ^ = p{C) , (3.29) 

namely 

(3.30) 



-1 



logp(C) 

Clearly if p(C) is very close to 1, i.e. if 

p(C) « 1 - 5 (3.31) 

with < S <C 1, we find 

5 1 - P(C) 



t~ ~ — ~ . . (3.32) 



One can also consider the inequality 

1-|A|<|1-A|, (3.33) 

which becomes an equality if the eigenvalue A is real and positive (since we are considering 
| A | < 1). This implies 

1 - n(C) = 1 - max |A| = min (1 - |A|) < min |1 - A| . (3.34) 

Xea(C) \€a(C) Xea(C) 

Thus, if p(C) is very close to 1 we obtain 

r « — — > — l — . (3.35) 

1 - p(C) mm Xea (c) |1 - A| 

Notice that eqs. (3.24), (3.26), (3.28) and (3.29) imply that, in the limit of large t, 

~ e-'/ r . (3.36) 



Thus, if we know the matrix C and we can find its eigenvalues A, we can evaluate the relaxation 
time r using eq. (3.30) [or the approximate expressions (3.32), (3.35)]. On the contrary, for a 
numerical determination of r one should use eq. (3.36) (for some definition of the norm 



4 Analysis of Critical Slowing-Down 

In this Section we review the analyses of critical slowing-down done in Refs. [8, 9] and we 
apply them to the four local algorithms considered in this paper. 6 The only difference with 
respect to those references is that, in our case, we minimize the free-field action (3.5) instead 
of thermalizing the configuration {f(x)}. Thus, their analyses can be applied directly to our 
case by setting the Gaussian noise rj (x) to zero. We stress that the results presented in this 
section are a straightforward application of the analyses reported in Refs. [8, 9] and that most of 

6 See also Ref. [15] for a careful analysis of CSD for the local hybrid Monte Carlo algorithm (which is equivalent 
to the stochastic overrelaxation method) applied to the free-field case using various updating schemes. Here we 
will consider only the so-called even/odd update. 
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these results, but not all of them, can be found there and in other articles. However, we believe 
that our presentation has some interesting insights, which can help the reader understand how 
local algorithms deal with the problem of critical slowing-down. Moreover, these results will be 
extensively used in Section 5, where we study generalizations of the overrelaxation algorithm and 
of the stochastic overrelaxation algorithm. 

Let us start by considering the overrelaxation update [see eqs. (3.16) and (3.18)] 

f° Ver \x) = - (w - 1 ) /(*) + £ E [ ]{x + e„) + f(x - e „) ] . (4.1) 

One can check that the condition {uj — l) 2 < 1 is sufficient to prove that this update never 
increases the value of the massless-free-field action [see eq. (3.22)]. Therefore we should have 7 
uo G (0,2). However, only when oo G (1,2) does one obtain [4] that the overrelaxation method 
performs better than the Los Alamos method, which corresponds to the case u — 1. For the 
Cornell method, it follows from eq. (3.23) that one should have aN(x) G (1,2), which gives 
a G (l/2d, 1/d) for the case (3 = oo. 

Clearly, using eq. (4.1), we need to know the value of the field / only at the site x and at the 
nearest-neighbor sites x ± in order to update f(x). Note that a site is defined to be even or 
odd according to whether the quantity 



d 



x 



E (4-2) 



is even or odd. Thus, if we consider lattices with an even number of sites in each direction and a 
checkerboard ordering, then we can first update all the even sites and subsequently all the odd 
ones, and so on. In order to implement the checkerboard update we can rewrite the update (4.1) 

as 



/mW = -(u-l)f t a (x) 

+ h Ed 1 + lft a (x + e») + mx-e,)] 
^ tt ft=i 

+ [fMx + eJ + f^x-ej]} , (4.3) 



where f a (x) are the three "color" components of f(x) and t denotes the number of sweeps through 
the entire lattice. Notice that when we update the odd sites, i.e. when ( — 1 )' x ' = —1, we use 
for the update the value of the newly updated field f° +1 at even sites. 
The idea in Neuberger's article [8] is to consider the Fourier transform 

f a (k) = E exp(-2mk ■ x) (4.4) 

X 

of eq. (4.3). To this end one can use the relation 

(-1) M = exp(-27riT • x) , (4.5) 

7 As shown in Section 3.3 of Ref. [4], this result also applies to the update given in eq. (2.3) when considering 
the minimizing functional (1.2). 
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where the vector T has components given by 

T 1 



/i — 1, . . . d 



■i* ~ 2 
In this way one obtains 

f~Uk) = -{u-l)j?{k) 

+ uc(k) [f?(k) + ft +1 {k) - f~?(k + T) + f? +1 (k + T) 

where we defined 

1 d 

c ( k ) = 2~rf ^ cos(2ttA; m ) . 



' fUk + T) j '{ f b a (k + T) 

where the 2x2 matrices A(k,u) and B(k,u) are given by 



-1 -1 
1 1 



Let us now define 



A(k,u) = 1_ + cjc(£;) 

= -{uj - 1) 1L + uc(k) ^ | ~J j 

C(ife,a;) = A- 1 (A;,cj)fi(A;,c«;) ; 



(4.6) 



(4.7) 
(4.8) 



One can write eq. (4.7) also for f? +1 (k + T), namely 
fUk + T) = -(cv -l)f~ t a (k + T) 

-ujc(k) [f?(k + T) + f? +1 (k + T) - f?(k) + f? +1 (k)} . (4.9) 
Then, equations (4.7) and (4.9) can be written as a system of two equations 



(4.10) 



(4.11) 
(4.12) 

(4.13) 

then, the update of the field f a , namely one even update followed by an odd update, can be 
written (in momentum space) as [see eq. (4.10)] 



f ~ Uk) )=c( M f* (fc) 

ft\i(k + T) ) \ f?(k + T) 



(4.14) 



Notice that the determinant of A(k,u) is equal to 1 for any u>, i.e. this matrix can always be 
inverted: 



A~\k,u) = i. + cuc(k) ^ _j _J j 



This gives 

C(k,u) = \2u 2 c 2 (k) -(w-l)ll + ujc{k) 



2 - u -u [1 + 2c(Jfe)] 
w [1 - 2c(ife)] -(2 - w) 



(4.15) 



(4.16) 
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It is easy to check that the eigenvalues of the matrix C(k,u) are given by 
\±(k,u) = 



2to 2 c 2 (k) - (w -l)]±yj [2uj 2 c 2 (k) - ( lu - 1 ) ] 2 - ( lu - 1 ) 2 (4.17) 

(4.18) 



(4.19) 



2uj 2 c 2 (k) - (lu - 1)J ± 2uu^Ju 2 c\k) - (lu - l)c 2 (k) 

Then, if we can prove that \\±(k,uj)\ < 1, we can use eq. (3.30) and write 

-1 



T 



log(max fc ^o \X±(k,uj)\) ' 



where we don't consider the constant (or zero) mode because it does not contribute to the action 
(3.5). 

We can obtain the eigenvalues \±(k,cu) also working in a slightly different way, namely 
following now Ref. [9]. The main difference with respect to the approach used above (based on 
Ref. [8]) is that, instead of considering the Fourier transform of eq. (4.3), one applies to eq. (4.1) 
the Fourier-like transform 



r*(k) = e r < 



X 



1 ± e 



-2-niT-x 



exp { — 2-nik ■ x) , 



(4.20) 



with T defined in eqs. (4.5) and (4.6). By using the result exp (±2iriT fl ) = —1 and the 
periodicity of the lattice one can verify that 



X 



-2niT-x 



1 ± e 

1 ± e ~ 2 * iT - x 



exp ( — 2nik ■ x 



exp ( — 2 7r i k ■ x) 



■ 2wik "f^(k) , (4.21) 



■ 2vik " f^(k) . (4.22) 



Thus, using eq. (4.1) and eqs. (4.20)-(4.22) and by updating first the f a ' + {k) components and 
then the f a '~(k) components, we obtain 



M{k,u) 



ft' + (k) 
ft'-(k) 



(4.23) 



M(k, w) = -(a; - l)l + 2w c(k) 



(4.24) 



where 

1 \ 

{uj - 1) 2uc{k) ) 

It is straightforward to check that the matrix M(k, uj) also has eigenvalues \±(k, uj) [given in eqs. 
(4.17) and (4.18)], namely M(k,cu) is the matrix C(k,u) [see eq. (4.16)] written in a different 
basis. Indeed, with 

(4.25) 



R = 



1 -1 



we have 



M(k,u) = RC(k,uj)R- L . 



(4.26) 



Note that we started the analysis of the overrelaxation update from one equation, namely 
eq. (4.1), and we are now considering a system of two equations. Thus, we should avoid double 
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counting for the momenta [9] and set, for example, kjN = 0,1, ...,N/2 — 1 and kiN = 
0,1, . . . , N — 1 for i — 1, . . . , d — 1. Therefore, in the limit of infinite lattice side N, the quantity 
c(fc) defined in eq. (4.8) takes values in the interval (—1/2, 1/2] and the lattice momentum p 2 (k) 
[see eq. (2.14)] takes values in [0, 4 d). Also, since we know that CSD is due to the long- wavelength 
modes, which are usually the ones with slowest relaxation, we should consider c(k) in the limit 
of small momenta. Obviously, if k^ = for all \i then c(k) = 1/2 and, using the fact that 
uo E (1,2) , one can check that \ + (0,uj) = 1 and A_(0,cj) = (1 — uo) 2 < 1, i.e. the constant (or 
zero) mode does not converge to zero. However, as said above, this mode does not contribute 
to the action (3.5) and it is therefore not relevant when studying CSD. On the contrary, for the 
smallest non-zero momentum — corresponding to ki — for % — 1, . . . , d — 1 and kd = 1/N - 
we have 

c sm (N) = ^[d- l + oos(^)] (4.27) 



and, in the limit of large lattice side N, 



1 - 



2tt 2 
dN 2 



= g [WW 



(4.28) 



namely we get the largest value of c(k) smaller than 1/2. This case is important because, as 
we will see below, the largest eigenvalue of the matrix C(k,uo) is very close to 1 exactly for the 
smallest non-zero momentum in the limit of large N. Thus, the quantity C(iV) will play a central 
role in the study of CSD for the four local algorithms considered here. 

Note that eqs. (4.27) and (4.28) are valid also in the case of asymmetric lattices if we set 
N = max^ N^. It is also interesting to check that c(k) switches sign when k^ goes to k^ + T M for 
all fj, [see eq. (4.8)]. Thus, we have | c(k) \ « 1/2 not only for small momenta p 2 (k) xs but also 
for the largest momenta p 2 (k) ps Ad, corresponding to small-wavelength modes. This is obvious 
if we observe that 



c{k) = 



1 



cos(2 7r/c M ) 



1 - 1 2 sin2 ( n K) 

a n=l 



1 




~ 2 


2d 



(4.29) 



where we used eq. (2.14). Indeed, if we consider the largest momentum p 2 (k), namely if we set 
ki — 1/2 for i — 1, . . . , d — 1 and kd — 1/2 — 1/N, then we obtain 



Clm(N) 



1 

2d 



— d + 1 + COS 71 — 



2tt 



1 

2d 



—d + 1 — cos 



2tt 



= -Csm(N) ■ 

(4.30) 

Therefore, in the limit of large lattice side N, the largest value of | c(k) \ smaller than 1/2 is given 
by 

\c(k(N)) \ = \\1- C(N)\ , (4.31) 

corresponding both to the smallest non-zero momentum and to the largest momentum. Since 
the eigenvalues X±(k, u) are functions only of c 2 {k) — or equivalently of | c(k) \ - and since they 
do not depend on the sign of the quantity c(k) [see eqs. (4.17) and (4.18)], the previous result 
implies that these large momenta contribute to CSD too. This unexpected effect is due to the 
even/odd update which couples the low- and high-frequency modes [15]. 
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4.1 Los Alamos Method 



As said before, the Los Alamos update is obtained from eq. (4.1) by setting ou — 1. Then, the 
eigenvalues X±(k,u) of the matrix C(k,u) become [see eqs. (4.8) and (4.17)] 



A_(£;,l) = (4.32) 

2 

X(k) = \+(k,l) = 4c 2 (fc) 



1 d 

- cos ( 2 77 M 



(4.33) 



and, in the limit of infinite lattice side N, we can consider kd taking values in the interval [0, 1/2) 
and ki in the interval [0, 1) for i — 1, . . . ,d — 1 [see comment after eq. (4.26)]. It is obvious that, 
if k^ = for all //, then A(fc) is equal to 1, while in any other case it is strictly smaller than I. 
Thus, all the Fourier modes relax, except for the constant (or zero) mode which, however, does 
not contribute to the action (3.5). Moreover, X(k) is always nonnegative. Therefore we can write 

< | \(k) | = X(k) < 1 . (4.34) 

It follows that [see eq. (4.19)] the relaxation time of the Los Alamos method is given by 

T = , , ~ l \77a\ (4-35) 

log(max fc7 , A(A;)) 

and if X(k) is very close to 1 — i.e. if \c(k)\ « 1/2 — we have [see eq. (3.32)] 

; (4.36) 



1 - max^o X(k) 

From the previous section we know that this is the case when one considers the smallest non-zero 
momentum — or the largest momentum — in the limit of large lattice side N. Then, using eq. 
(4.31) we obtain 

max X(k) = [1 - ((N)] 2 » 1 - 2((N) , (4.37) 

k =£ 

from which follows 

TLosAlarnos ~ J^J]^ = ^ " ( 4 - 38 ) 

So, as expected, the dynamic critical exponent z is equal to 2. 
4.2 Overrelaxation Method 8 

As seen above, the update for the overrelaxation algorithm is given by eq. (4.1) with the parameter 
ui taking values in the interval (1, 2) and its relaxation time is related to the eigenvalues X±(k, u) 
of the matrix C(k,u) [see eqs. (4.17) and (4.18)] through the relation [see eq. (4.19)] 

~ 1 (4.39) 



log(max fc7 , |A±(fc)|) ' 



The results presented in this section apply also to the Cornell method [see eqs. (3.17), (3.18) and (3.23)]. 
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provided that max^ |A±(/c)| < 1. As said at the end of Section 3, we must verify that this 
condition is satisfied to ensure the convergence of the algorithm. To this end, following Ref. [8], 
we define 



_ 2uj 2 c 2 (k) 
uj - 1 

and, using the fact that uj G (1, 2), we can write eq. (4.17) as 



rM = 1 _ — - ^ (4.40) 

UJ — 1 



\±(k, uj) — (uj — 1 ) 



— r(k,u) ± \J r 2 (k, uj) — 1 



(4.41) 



Note that r(k,u) G (— oo, 1]. There are therefore two possibilities, described below. 



1) r(k,uj) G [—1, 1] for all values of k, namely we have to impose the condition r(k,uo) > — 1 
for all k. In this case the eigenvalues \±(k,cu) are complex conjugates of each other and 
we have 

| \±(k,oj) | = \u - 1| , (4.42) 
namely |A± is independent of k (see Refs. [2, 8, 9]). Also note that 

\u-\\=u-\<\. (4.43) 

2) r(k, uj) < —1 for some values of k, denoted as k. In this case the corresponding eigenvalues 
\±(k,u) are real and we can verify that < \^(k,uj) < X + (k,uj). From eq. (4.40) it is 
clear that we can obtain r(k,uj) < —1 only for the largest values of c 2 (k), i.e. c 2 (k) 1/4, 
corresponding to momenta p 2 (k) ~ or to momenta p 2 (k) rs Ad (see comment at the end 
of Section 4). Considering now eq. (4.40) with r{k,uj) < —1 we get 

< uj - 1 < uj 2 c 2 (k) (4.44) 

and 

uj 2 c 2 (k) < 2uo 2 c 2 (k) - (uj - 1) . (4.45) 
Also, since c(k) G (—1/2, 1/2] for all values of k, we have 

c 2 (k) < \ (4.46) 



and one can check that 



and 



< 2uj 2 c 2 (k) - (uj - 1) < Uuj - I) 2 + \ (4.47) 



uj 2 c 2 (k) - (uj - 1) < (2 A U)2 . (4.48) 



These inequalities imply [see eq. (4.18)] 

\ + (k,uj) < l -(u - I) 2 + l - + uj^^- < 1 ; (4.49) 

moreover, X + (k,uj) is equal to 1 only if c 2 (k) = 1/4, i.e. if k = 0. Finally, from eqs. (4.17), 
(4.44) and (4.45) one can prove that \ + (k,uj) > uj — 1. 
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To sum up we can say that, with u G (1, 2) , there are two different situations. If the inequality 
r(k,uj) > —1 is satisfied for all values of k, then the eigenvalues \±(k,u) are complex conjugates 
of each other with absolute value independent of k and given by uo — 1 . In this case we have 
that the largest eigenvalue is close to 1 — which is the interesting case when one studies CSD 
- only when uo is close to 2. If, on the contrary, we have that r(k,uo) < — 1 for some values k, 
then the largest eigenvalue of the matrix C(k,uo) is real and given by max^ \+(k, ui). (Note 
that max^Q is taken only over the values k.) In this case it is easy to check that, in the limit of 
large N, 

max A + (JU) - 1 - = 1 - - r2 , (4.50) 

Mo 2 - uo (2 - uo)dN 2 ' v ; 

where we used eqs. (4.17) and (4.31). 

In order to study CSD for the overrelaxation algorithm one usually writes 

u = rh (4 - 51) 

with fl G (0, 1). By using eq. (4.51) we can then rewrite r(k,cu) as 

r(h,Q) = 1 - . (4.52) 

We now discuss separately the two cases considered above. 

Case 1), i.e. if r(k,cu) > —1 for all values of k, corresponds to 

Q 2 < 1 - 4c 2 {k) . (4.53) 

Since this condition should be satisfied for all values of k we must have 

VI 2 < 1 - 4 max c 2 (k) , (4.54) 

where again we don't consider the zero mode k — 0. When N goes to infinity we find [using eq. 
(4.31)] 

^ 2 <2C(A0 = ^, (4.55) 
namely fl goes to zero. In this case we have 

|A±(fc,o;)| = uj-1 = (4.56) 
In the limit of small Q we obtain 

| \±(k,w) | w 1 - 2Q (4.57) 

and the relaxation time becomes [see eq. (3.32)] 

1 1 



max fc ^ | \±(k,u) \ 2Vt 



(4.58) 
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So, if we want to minimize the relaxation time r we have to maximize the value of f2 allowed by 
the inequality (4.55), i.e we have to set 

n = -4^~ • (4-59) 
VdN 

In this way we get 

Tover ~ ^ N (4.60) 
4 7T 

and, as expected, we have z — 1. Note that the tuning for Q given by eq. (4.59) implies 

u = = 2 (l + 1 . (4.61) 

Case 2) corresponds to the existence of values k satisfying the inequality 

n 2 > 1 - Ac 2 (k) . (4.62) 

In this case, as we saw above, we have 

maxA + (M~l- (2 _ 4 ^ 2 - (4-63) 



This implies 

1 (2-u)dN 2 



(4.64) 



1 — max^Q \ + (k,oj) Aujtt 2 

and in order to minimize the relaxation time r we have to minimize (2 — uo)/uo = Q. This can 
be done by choosing Q ~ 1/N m , yielding r proportional to N 2 ~ m and z = 2 — m. In particular, 
it might seem possible to set m = 2 so that r becomes constant in iV and CSD is completely 
eliminated. However, we note that the tuning of Q must be done while still satisfying eq. (4.62), 
since it defines what we are calling "case 2)". This implies the condition 

n 2 > l-4c 2 (k) ~ 4?, (4.65) 



_/V 2m -- - V , jy 2 

namely m < 1 and z > 1. Thus, since we want to minimize Q, we have to set m — 1 and impose 
that the inequality (4.62) become the equality 

Vi 2 = 1 - 4 max c 2 (k) . (4.66) 

We can conclude by saying that — in both cases considered above — the best tuning for 
the overrelaxation algorithm is obtained from the condition (4.66). Then, the largest eigenvalues 
X±(k,u) of the matrix C(k,u) with k ^ are real and both equal to u — 1 [corresponding to 
r(k, Q) = —1], while for all the other non-zero momenta these eigenvalues are complex conjugates 
of each other [corresponding to r(k, Q) > —1]. 

If we do not tune the value of the parameter Q (or equivalently of the parameter a;) , we can 
consider two limiting cases. When is so small that the condition 

VI 2 < 2((N) (4.67) 
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is satisfied for all the lattice sides N considered, then r(k,cu) > — 1 for all momenta and we 
obtain [see eq. (4.58)] 

r » ^ • (4.68) 

Thus, the relaxation time r is constant in N and we get z — 0. However, r is very large, and 
even though the number of gauge-fixing sweeps n g j is in this case independent of the lattice 
side N, we need a very large n g f in order to complete the gauge fixing even for relatively small 
lattice volumes. On the contrary if Q is large and, for all lattice sides N considered, we can find 
non-zero momenta p 2 (k) such that the condition 

VI 2 > 1 - Ac 2 (k) (4.69) 

is satisfied, namely r(k,u>) < —1, then the relaxation time r is given by [see eq. (4.64)] 

(2 - u)dN 2 



4o> 7T 



2 



(4.70) 



and z = 2. Nevertheless, in this case, the overrelaxation algorithm works better than the Los 
Alamos method: in fact, using eq. (4.38), we obtain 

T ~ T~LosAlamos = ^TLosAlamos (4-71) 

and r is always smaller [for u G (1,2), i.e. f2 e (0,1)] than the Los Alamos relaxation time. 
Clearly, when N goes to infinity with Q fixed, we always obtain eq. (4.70) and find z — 2. 

As a final remark, we note that this analysis implies that the dynamic critical exponent z of 
the overrelaxation algorithm depends only on the relation between Q 2 and 2((N), namely if for 
all values of N considered we set Q 2 much smaller than, equal to, or much larger than 2 C(iV) we 
have that the dynamic critical exponent z is equal to 0, 1 or 2 respectively. 



4.3 Stochastic Overrelaxation Method 

In this section we want to analyze the critical slowing-down of the stochastic overrelaxation 
method. As explained in Ref. [4], this algorithm is similar in spirit to the idea behind the so- 
called hybrid version of overrelaxed algorithms (HOR), which are used to speed up Monte Carlo 
simulations of spin models, lattice gauge theory, etc. [2, 9, 16]. In these algorithms, m micro- 
canonical (or energy conserving) update sweeps are done followed by one standard local ergodic 
update (such as Metropolis or heat-bath) sweep over the lattice. Actually, for the Gaussian 
model, it has been proven [9] that the best result is obtained when the micro-canonical steps 
are chosen at random, namely when m is the average number of micro-canonical sweeps between 
two subsequent ergodic updates. This is essentially what is done in the stochastic overrelaxation 
method [see eq. (2.4)], with mj [m + 1) equal in average to p or, equivalently, m equal on average 
to pj (1 — p). 

In order to follow the analysis in Ref. [9], we suppose that the stochastic overrelaxation 
method is implemented as an HOR algorithm: m sweeps using the "micro-canonical" update 
[see eq. (4.1) with u — 2] 

J(x) = - J(x) + jE[/> + ^) + 7(* - e„) ] , (4.72) 

a n=l 
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which does not change the value of the minimizing functional, and one sweep using the Los 
Alamos update [see eq. (4.1) with u — 1] 

/(*) = E [ 7(* + + 7(* - e .) ] . ( 4 - 73 ) 

which brings the minimizing functional to its local absolute minimum. Let us notice that, for 
u; — 2, the eigenvalues \±(k,u) in eq. (4.18) become 



A±(Jfe,2) = [-1 + 8c 2 (/c)J ± i4^c 2 (k) - 4c 4 (/c) = exp[±i9(k)] , (4.74) 
where we define 6{k) such that 

cos [#(£;)/ 2] = 2c(fc) . (4.75) 

Clearly, in this case, we have |A±(fc, 2)| = 1 for all values of k and, as expected, none of the 
Fourier modes relaxes. 

If we consider as one sweep of the lattice the combination of m sweeps using the update (4.72) 
and one sweep using the update (4.73), then the matrix that defines this combined update (in 
momentum space) is given by 

C(k,m) = C(£;,l) [C(k,2)] m , (4.76) 

where C(k,cu) is defined in eq. (4.16), or by 

M(k,m) = M(k,l) [M{k,2)} m , (4.77) 

with M(k,u>) defined in eq. (4.24). Since it is easier to work with M(k,u) than with C(k,u>) we 
will use eq. (4.77). However, one can check that C(k, m) and M(k, m) have the same eigenvalues 9 . 
Then, following Ref. [9], we can write 

where 6{k) is defined in eq. (4.75) and the matrix V(k) is given by 

„ (i)s /«p[i«(t)/2] «p[-«(t)/2]\ _ (480) 

This implies that 

, / exp[-im0(ifc)] \ w x . . 



9 This is immediate if we consider that the matrix R defined in eq. (4.25) is independent of k and u>. Therefore 
the relation (4.26) between the matrices M(k,u)) and C(k,w) is also valid for the matrices M(k,m) and C{k,m). 
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One can also write the matrix M(k, 1) in the dyadic form 

' 1 



M(k,l) = 2c(k) 



2c(k) 

and, using eqs. (4.77), (4.82) and (4.83) we obtain 
2c(k) ( 1 \ 



= 2c(k) 



2c(k) 



(0,1) 



(4.83) 



M(k,m) = — 



sm[6(k)/2] \ 2c(k) J 
The eigenvalues of this matrix are equal to zero and to 



sin [m 6(h) ], sin 



m+-)9(k) 



(4.84) 



cos 



0(k) 



COS 



m+ - ) 9(k) 



(4.85) 



However, m is not fixed but varies between zero and infinity with probability p m (1 — p). 
This gives an average value 



P 



(m) = E m P m (l -P) = h*g I = pi log Ef^ — 

2^m=0 P a P m =0 1 P 



(4.86) 



m=0 



as said above. Thus, instead of looking for the eigenvalues of the matrix M(k,m) we should 
consider the matrix 

oo 

(4.87) 



M(k,p) = (1 - p) p m M(k,m) 



m=0 



After writing eq. (4.84) as 
2c(k) 



M(k,m) 



x 



2c(k) 



Im ( — exp[im6(k)]i exp 



% [m+ - ) 9(h) 



, (4- 



it is straightforward to check that the matrix Ai(k,p) can be written in the dyadic form 

(1 - p) cos [#(£;) /2] 



M(k,p) = 



[1 + p 2 - 2p cos9(k)} sm[9(k)/2' 
( 1 



x 



cos [#(£)/ 2; 



(-p sin9(k), (1 + p) sin[9(k)/2\) (4.89) 



and has eigenvalues zero and 

A(M 



(1 - pf cos 2 [9(k)/2] 
(1 - p) 2 + 4psin 2 [^(A;)/2; 



(4.90) 



Clearly this eigenvalue is nonnegative for any p e (0, 1) and for p = (i.e. ( m ) = 0) we obtain 
the non-zero eigenvalue A(/c,0) = cos 2 [#(A;)/2] = 4c 2 (/c) of the Los Alamos method [see eq. 
(4.33)]. 
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Note that, since the matrix M.(k,p) describes (on average) (to) + 1 sweeps of the lattice, 
eq. (3.29) is not correct for the stochastic overrelaxation algorithm. In fact, in this case, the 
relaxation time r is related to the eigenvalue X(k,p) by the expression [9] 

max | X(k,p) | = max X(k,p) = e - (<m) + 1)/T , (4.91) 

namely 

T -1 

! = : . (4.92) 

(to) + 1 log(max fe ^o X(k,p)) ' 

In order to study CSD for the stochastic overrelaxation algorithm, we introduce P G (0,1) 
and write 

P = \^ ■ ( 4 - 93 ) 
Then, the eigenvalue in eq. (4.90) becomes 

M k p) = p2 cos W)/2] (4 94) 

Al/C '^ j 1 + (P 2 - l)cos 2 [W/2] • l4 ' y4j 
Also, using eqs. (4.29) and (4.75), one obtains 

cos[0(fc)/2] = 2c(k) = 1 - ^ ee 1 - f d , (4.95) 

where r is the magnitude | p(k) | of the lattice momentum. Note that, since c(k) takes values in 
the interval (-1/2,1/2], we have cos[6(k) / 2} e (-1,1], 0(k) E [0,2 n) and r G [0,2y/d). Thus, 
we can rewrite this eigenvalue as 

A(r, P) = ^ 2d) . (4.96) 

1 + (P 2 - 1) (l - fi) 

It is easy to check that 

• A(r,P) = 1 for r = 0, 

• the derivative of A(r, P) with respect to r is zero for r = and r = \^2d, 

• this derivative is negative (respectively positive) for r < y/2d (respectively for r > V2d). 

In Fig. 1 we plot A(r,P) as a function of r for the case d — 4 and P = 0.2. We note that 
the eigenvalue A(r, P) does not show the oscillatory behavior that can be observed when one 
considers a probability distribution that is uniform 10 for to in the interval [1,2m -1] (see Figs. 
1 and 2 in Ref. [9]). Also note that this eigenvalue is close to 1 for small momenta r and 
for very large momenta, i.e. with r xs 2\fd ) and that in both cases we have c 2 (k) w 1/4 and 
cos 2 [0(fc) / 2] f« 1 . As explained at the end of Section 4, this is a natural consequence of the 
even/odd updating scheme. 
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Note that with this distribution we automatically have ( m ) equal to m . 
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It follows that maxfc^o A(/c, P) < 1 is obtained for c 2 (k) ps 1/4. Then, using eqs. (4.31) and 
(4.75) we have, in the limit of large lattice side N, 



cos 2 [8(k)f 2] = Ac 2 (k) « 1 - 2C(iV) 



and from eq. (4.94) we find 

A(fc,P) w 1 
Thus, X(k, P) is very close to 1 and we obtain 

T 



2C(AQ 
P 2 



(to) + 1 1 - max fc ^o A(fc, P) 
Using again eq. (4.94) and eq. (4.97) we get 

2C(AQ 



l-maxA(A;,P)^^ riv)(i _ p2) + p2 



so that 11 



(to) + 

Also, from eqs. (4.86) and (4.93), we have 



T ^(l-P 2 ) + 



2((N) 



(to) 



V 



1 -P 



and 



These equations give 



(to) + 1 



(l - P 2 ) + 



1-p 2P 
1 + P 



2P 

P 2 
2C(iV) 



1 + P 
2P 



We can now fix P by minimizing the value of r. In this way we obtain 

2C(A0 ^ 2C(iV) 



(4.97) 
(4.98) 

(4.99) 



(4.100) 
(4.101) 

(4.102) 

(4.103) 
(4.104) 

(4.105) 



[1 - 2C(A0] (1 + 2P) (1 + 2P)' 
Therefore, in the limit of large lattice side N, we get that P goes to zero (and p goes to 1) as 

2tt 



and 



Tstoc 



1 + 



P 2 



2((N) 



1 1 VdN 



2P 



2tt 



(4.106) 



(4.107) 



Since at this point we don't know the relation between P and Q{N) we have to keep all terms in this equation. 
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Thus, as expected, we have z — 1. Let us notice that the tuning for P given in eq. (4.106) 
coincides with the tuning for Q obtained in eq. (4.59) for the overrelaxation algorithm. Also we 
have 

1 - P 47T 

p = « 1 - 2P = 1 =— (4.108) 

and 

1 + " = TTp " TTH = "' (4109) 

namely we find the relation p m — This is in agreement with the result obtained analytically 
and numerically in two dimensions at finite f3 (see Sections 5 and 7.3 in Ref. [4]). 

Finally, if we do not tune the value of the parameter P (or equivalently of the parameter p) 
then we have again two limiting cases. In fact, if P is very small and such that 

P 2 < 2((N) (4.110) 

then from eq. (4.104) we obtain 



1 

2P ' 

namely the relaxation time is constant in N and z — 0. On the contrary, if P is large and 



(4.111) 



P 2 > 2((N) , (4.112) 

then we find 

1 + P P 2 _ (l + P)P dN 2 
T ~^2P~2C(N)- 2 4^ (4113) 

and [using eq. (4.38)] 

(1 + P) P 

t ~ T LosA i amos . (4.114) 

This gives a dynamic critical exponent z equal to 2. However, in this case, the improved local 
algorithm without tuning works better than the Los Alamos method: the relaxation time r given 
in eq. (4.114) is always smaller [for P e (0, 1)] than the Los Alamos relaxation time. As in the 
overrelaxation case, when goes to infinity (with P fixed) we always get (4.113) and find z — 2. 
Also, this analysis implies that the dynamic critical exponent z depends only on the relation 
between P 2 and 2£(AT). In fact, if we set P 2 much smaller than, equal to, or much larger than 
2((N) we have that the dynamic critical exponent z is equal to 0, 1 or 2 respectively. 



5 Generalized Local Algorithms 

In Sections 2, 3 and 4 we have studied CSD for two main types of local algorithms: the over- 
relaxation algorithm — which coincides with the Los Alamos algorithm for uj — 1 and with the 
Cornell algorithm with the choice [see eq. (3.23)] uj ps aj\f(y) — and the stochastic overrelax- 
ation algorithm. In both cases we have seen that with a careful tuning (of the parameters uj and 
p respectively) we obtain a dynamic critical exponent z — 1. In this section we want to see if it 
possible to generalize these algorithms in order to get z < 1. To the best of our knowledge, the 
results presented in Sections 5.2 and 5.3 below are new. 
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5.1 Generalized Overrelaxation Algorithm 

Let us consider the following generalization of the standard overrelaxation update [see eq. (2.3)] 

g(Hnear)^ = d(uj)h)(x) + bjuj) g(x) 

y/[a(u) + b{u) f - a(u) b{u) [2 - T{x) ] ' 

One can check that g( hnear )( x ) e SU{2) and that the overrelaxation update corresponds to the 
choices a(u) = u and b{oj) = 1 — u. In the limit of large number of gauge-fixing sweeps t [see 
Section 3] we can use eq. (3.12) and obtain 



g( linear \x) 



cl(uj) ~, b(uo) 



I a(w) + 6(w) | | a(w) + b(u) | w v ' 

where terms of order e 2 have been neglected. This udate can be written in the simpler form 12 
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(linear) 



x) ?s a{uo) h\x) + b(u)g(x) (5.3) 



if we assume the condition a (a;) + b(u) = 1 or, equivalently, by a redefinition of the coefficients 

— * 

a{uj) and b(u). Finally, using eq. (3.4), we find the update for the f(x) field 

f mear \x)=a(.)f L ° sAL \x) + b(u J )M. (5-4) 

— * 

If we consider the dependence of the massless-free-field action on the value of the field / at a given 
site y [see eqs. (3.21) and (3.22)], it is clear that this is the most general local linear update (with 
site-independent coefficients) of the field f(y). One can also check that the condition b 2 {uj) < 1 
is sufficient to prove that the update (5.4) never increases the value of the massless-free-field 
action [see eq. (3.22)]. Since the definition of uj is in principle arbitrary, we can at this point 
set a(u) = u and obtain that the standard overrelaxation algorithm is the most general local 
linear update with side-independent coefficients. Therefore, following the analysis presented in 
the previous section we have at best that z — 1. 

In order to understand why, in this case, one cannot get a dynamic critical exponent z smaller 
than 1, we can follow Ref. [10] and consider the inequality [see eq. (3.35)] 

T ^ rr~~ \ IT- m • ( 5 - 5 ) 

mm^o I 1 - A±{k,u>) \ 

When r 2 (k,u) < 1, namely when these eigenvalues are complex, one finds [see eqs. (4.40) and 
(4.41)] 



1 1 - \±(k, u)\ = uj^Jl - A c 2 (k) . (5.6) 
Since the tuning condition, i.e. r 2 (k,u) < 1, is equivalent to the relation [see eq. (4.53)] 

Q < yfl - 4c 2 (k) (5.7) 

[where the inequality becomes an equality for the largest value of c 2 (k) < 1/4], we have [using 
eq. (4.51)] 

mjn|l- A ± (*,a;)| = ^ • (5-8) 



12 



We have considered this type of update also in Section 5 of Ref. [4] . 



25 



In the limit of large N this implies [see eq. (4.31)] 



A < y/2C(N) (5.9) 

and 

mjn|l-A ± (*,a;)| « 2^2((N) = -ijL , (5.10) 



so that [in agreement with eq. (4.60)] 



(Ml) 



which gives z — 1. 

Equations (5.6) and (5.10) are the starting point of Ref. [10] (see their Fig. 1). The idea 
in that article is that one should look for an update characterized by a matrix C(k) whose 
eigenvalues X(k) satisfy, in the limit of large N, the relation 

mm | 1 - X(k) | oc ^_ j (5.12) 
with m > 2. In fact, this would imply 

r ~ N 2/m , (5.13) 

namely z < 1. From eq. (4.41) it is clear that we have m = 2, and therefore z — \, because the 
eigenvalues \±(k,u) are solutions of a quadratic equation, i.e. because C(k,u) is a 2 x 2 matrix. 
What is proven in Ref. [10] is that, unfortunately, we cannot get m > 2 because any general 
update matrix C can be shown to be block-diagonal with blocks of size not larger than 2x2. In 
the next section we will consider explicitly local algorithms characterized by updating matrices 
of size 4x4 and we will check that we cannot obtain z smaller than 1 in that case. 



5.2 A 4 x 4 Updating Matrix 

Following what was done in Section 4, we consider here the update (4.1) and the Fourier-like 
transformation defined in eq. (4.20) and generalize that procedure in a way that produces an 
updating matrix that is 4 x 4 instead of 2 x 2. To this end, let us recall that if the vector T has 
components T M = 1/2 for all \x [see eq. (4.6)] then we have 

, . ._ x f 2 if I x I is even , , N 

1 + exp(-2mT ■ x) = I \ \ . 5.14 

v ' [ if | x | is odd v ' 

1 — exp ( — 2-kiT ■ x) = < ^ !p!" r !! S ev i e * 1 (5.15) 
v ' [2 if | x | is odd v ' 

These two linear combinations automatically select the even and odd sub-lattices and can be 
used to construct the two-component field f b,± (k) [see eq. (4.20)]. Let us notice that one can 
also rewrite the above relations as 

exp ( — 2-aiTx ■ x) ± exp ( — 2-niT 2 • x) (5.16) 
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with Ti )A( = 1 and T 2iM = 1/2 for all //. We can easily generalize this result and divide the lattice 
in, for example, four sub- lattices. 13 In fact with T 1:/J/ = 1, T 2 ^ = 1/2, T 3i ^ = 3/4 and T 4iM = 1/4 
we have that the linear combinations 

L ee ( X ) = e^-^iT,-*) + e (-2*iT 2 . X ) + e (-2mT 3 - X ) + e (-2 ff iT 4 .x) (5^7) 

L eo (x) = e ( - 2 ^ Tl '^ + e (- 2 ^-x) _ e (-2 7 r J T 3 .x) _ e (-2^T 4 .x) ^g) 

L oe ( X ) = e^^m-x) _ e (-2viT 2 . X ) + ie (-2niT 3 . X ) _ ie (-2niT 4 . X ) (5JQ) 

L 00 (x) = ei-^^.x) _ e (-2*iT 2 . X ) _ ie (-2niT 3 . X ) + ie (-2niT 4 . X ) ^ ^ 2 0) 

called respectively even-even, even-odd, odd-even and odd-odd, are always zero but for the 
following cases 

L ee (x) = 4 if |x|mod4 = (5.21) 

L eo (x) = 4 if \x\ mod 4 = 2 (5.22) 

L oe (x) = 4 if |x|mod4 = l (5.23) 

L 0O (x) = 4 if | x | mod 4 = 3. (5.24) 

Thus, in this way we can automatically divide the lattice in four sub-lattices and define a four- 
component field /-(fc), /•(*), /• (k)) by 

= E f a (x)Lee(x) exp(-2mk ■ x) , (5.25) 

and analogously for the other three components. Also, in analogy with eqs. (4.21) and (4.22), 
one can check that 

E Fix + e^L^x) exp(-27rik ■ x) = e +2 * ik » /»(£;) (5.26) 

X 

E f a (x-e f ,)L ee (x)exp(-27rik ■ x) = e" 2 ^ f a oe (k) (5.27) 

X 

and similarly using the other linear combinations defined in eqs. (5.18)-(5.20). Then, for any 
updating sequence of the four components (fe e (k),f^ (k),f^ e (k),f^ (k)) we obtain an updating 
matrix which is 4 x 4. For example, if we update these components in the order even-even, 
even-odd, odd-even and odd-odd one can verify that the updating matrix is given by 14 

( E~(k) \ 

(1 - w)L + w U-^ + m .w-wJ' (5 ' 28) 

where 1L is the 4x4 identity matrix and E ± (k) are 2x2 matrices defined by 



13 In this section, in order to simplify the notation, we suppose that the lattice side N is a multiple of 4. 

14 In order to avoid over-counting for the momenta we should set, for example, kd N = 0,1, . . . , iV/4 — 1 and 



hN = 0,1, ...,N- 1 for i = 1, ...,d- 1. 
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with 



,±2 7ri fc„ 



= e r (fc) ± i ei(k) . 



(5.30) 



Clearly the matrix in eq. (5.28) has a structure very similar to that of the 2x2 matrix given in 
eq. (4.24). In fact, in this case, we are still using a checkerboard ordering (first all even sites and 
then all odd sites), but with the even and odd sub-lattices divided in turn into two sub-lattices. 
The eigenvalues of this 4x4 matrix are, as expected from Ref. [10], solutions of two different 
second-order equations and are given by 



\ r ± (k,uu) = (1 - uj) + 2uj 2 e 2 {k) ±2uyJ{l- uj) 2 e 2 {k) + uj 2 e A r (k) (5.31) 
Xi{k,uj) = (1 - uj) + 2uj 2 e 2 {k) ± 2uJ{l - uo) 2 e 2 {k) + uj 2 e\(k) . (5.32) 



Since e r (k) is equal to the quantity c(k) defined in eq. (4.8), we have that the eigenvalues X±(k, uj) 
coincide with the two eigenvalues of the overrelaxation method obtained in Section 4. Also, for 
uj — 1, we have \±(k,uj) = 0,4e^(/c) and 4e 2 (k) = 4c 2 (/c) is the non-zero eigenvalue of the Los 
Alamos method found in Section 4.1. We note that 



= 2d ^ sin ( 27T K) 



(5.33) 



assumes its largest value (equal to 1/2) when k^ = 1/4 (for all directions /i). This implies [see 
eqs. (5.31) and (5.32)] that CSD is now due not only to the long- wavelength modes — i.e. small 
momenta p 2 (k) w 0, for which 4c 2 (k) ~ 1 — but also to the modes with k^ pa 1/4, corresponding 
to momenta p 2 (k) pa 2d and for which 4e 2 (k) ~ 1. This result can be explained by observing 
that the division of the lattice in four sub-lattices couples the modes with k^ pa to the modes 
with k^ pa 1/4. In fact, e^k) becomes e r (k) = c(k) when k^ goes to k^ + 1/4 (for all directions 

/')• 

Of course, if one considers a different updating sequence for the four components of the field 
f a (k), then the updating matrix will be different from that reported in eq. (5.28). For example, 
if we update these components in the order even-even, odd-even, even-odd and odd-odd one can 
verify that the updating matrix is given by 



M 4 



where 



/ (1 - uj) uje~(k) uje + (k) \ 

(1 - uo)uo 2 e +2 {k) f(k,uj) uoe+(k) f(k,uj) g{k,uj) 

(1 - uj)uje+{k) we-(fe) f(k,uj) uj 2 e +2 (k) 

\ (1 — Lo)g(k,uj) uoe + (k) f(k, uj) h(k,uj) f(k,uj) + uj e + (k) g(k,u>) J 



f(k,uj) = (1 - uj) + uj 2 e + (k)e-(k) 
g(k,uj) = uje~(k) +uj 3 e +3 (k) 
h{k,uj) = uj 2 \e~ 2 {k) + e + \k) f{k,uj) 



, (5.34) 



(5.35) 
(5.36) 
(5.37) 
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In this case we were not able to prove that the characteristic equation of the above 4x4 matrix can 
be factorized in two different second-order equations and we could not find any simple expression 
for the four eigenvalues. [On the other hand, from (5.34) it is obvious that one of the eigenvalues 
of the matrix Ai^ is zero when ou — 1.] However, after some manipulations one can verify that 
the characteristic equation of can be factorized into two "almost" second-order equations 
given by 



. 



(1-lu-\) 2 - 2uj 2 \e + (k)e~(k) ± V\u 2 [e (k) + Ae + (k) 
By using eq. (5.30) and the relation e r {k) = c(k) one can re-write the above equations as 
{1-uj-Xf - Auj 2 Xc 2 {k) + V\u 2 \c 2 {k) (l + V\) 2 

-e 2 {k) (l + V\) 2 - 2ic(k) ei (k)(l - A) 



(5.38) 



(5.39) 



and 



[1-uj-X) 2 - Auj 2 Xc 2 (k) - V\u 2 \c 2 (k) (l - V\) 2 

-ej(k) (l - V\f - 2ic(k) ei (k)(l - A) 



. 



(5.40) 



Let us notice that, written in this way, these two equations are very similar to the characteristic 
equation of the matrix M(k,u) given in eq. (4.24) [or equivalently of the matrix C(k,u) given 
in eq. (4.16)]: 

(1 - u - A) 2 - Au 2 \c 2 (k) = . (5.41) 

Before considering the equations (5.39) and (5.40) it is interesting to study in more detail eq. 
(5.41) and see how we can estimate the dynamic critical exponent z. To this end we can re-write 
the previous equation as 



+ 2uj (1 - A) { Jl - 4c 2 (k) - 1 - 2uc 2 (k) \ = (5.42) 



(1 - A) - ujjl - 4c 2 {k) 

and for the largest value of c 2 {k) < 1/4, in the limit of large lattice side N, we obtain [using eq. 
(4.31)] 



\l-\)-uj20N) + O(C 3 / 2 (N)) 

+ 2u(l- A) {v/2C(A0 - (l - |) + 0(C(N))} = . 
Then by setting A = 1 — 5 we get 

5 2 + 2u 2 ((N) + 2^(|-l) + 0(5((N)X 2 (N)) = . 
Thus, if we do not tune the parameter uj we have that 5 is the solution of the equation 

u((N) 0(5 2 ,5aN)X 2 (N)) = , 



(5.43) 



(5.44) 



(5.45) 
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namely 5 oc C(^0 ~ l/N 2 . This implies [see end of Section 3] 

r > 7^7 ~ N 2 (5.46) 
I ° I 

and z — 2. On the contrary if we tune cu — 2/(1 + fi) by setting 

n = n c m (^) , (5.47) 

with m > 0, then from eq. (5.44) we obtain 

5 2 - 45QC n {N) + 8£{N) + 0{5 ({N),5 ( 2m {N),( 1+m {N)),( 2 {N)) = . (5.48) 

Note that for m < 1/2 the previous equation simplifies to 

5 2 - 45QC m (N) + 0(5( 2m {N),({N)) = , (5.49) 

with solutions S = and 8 = 4Q( m (N). Thus, in this case one of the eigenvalues is equal to 1 
and the algorithm does not converge. On the other hand, for m = 1/2 we have 

5 2 - + 8((N) + 0(5C(N),C 3/2 (N)) = (5.50) 



and the solutions 5 are clearly proportional to y£(iV), giving z — 1. One can also check that 

the tuning condition f2 = y/2 [see eq. (4.55)] reduces the previous equation to a perfect square, 
namely 

2 + 0(5C(N),C 3/2 (N)) = 0. (5.51) 



5 - 2pC(N) 

Clearly, this analysis can also be applied to the two equations (5.39) and (5.40). In this 
way we will check that the dynamic critical exponent z cannot be smaller than 1 when using 
the updating matrix Ma. To this end, we have to consider the largest value of c 2 (k) < 1/4 - 
corresponding to ki — for i — 1, . . . , d — 1 and kd — l/N — in the limit of large lattice side N, 
i.e. we use eq. (4.31). For the same k, in the limit of large N, we also have the relation 



^ - hi Sm &) » ^ + - ^ + 0(C 3 ' 2 (N)) . (5.52) 

Then, for = oo and by considering the equation (5.39) we get 

(1 - A) 2 + 2cu (1 - A) (| - l) + V\ ^ (l + V\) 2 = . (5.53) 

Obviously, with uo G (1,2) there is no solution A ~ 1 that can satisfy this equation. In other 
words, if we set A = 1 — 5, then 5 stays finite and we don't get a divergent relaxation time r. 
Thus, critical slowing-down should be related to the solutions of the equation (5.40). For this 
equation, using eqs. (4.31) and (5.52) and setting again A = 1 — 5, we find 



(l-^) 5 2 + 2c 2 C(A0 + 2c^(!-l) 



+ tu 2 5 J C -^T + 0(5\5 2 JaN),5C(N)X 2 (N)) = . (5.54) 



2d 
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This is very similar to eq. (5.44) and we can check that the two new terms — uj 2 5 2 /1Q and 
i J 1 5 ^J((N)/(2d) do not spoil the analysis already done for that equation. In fact, if we do not 
tune the parameter uj we have 

u((N) +6(^-1)+ 0(5 2 ,5j((N),( 2 (N)) = , (5.55) 

namely 8 oc ({N) ~ 1/N 2 and z — 2. On the contrary, if we use the tuning given in eq. (5.47) 
and consider < m < 1/2 we have 

^5 2 - 45U( m {N) + 0{5 3 ,5 2 ( m {N),5( 2m {N),({N)) = (5.56) 

with solutions 5 = and 5 = 16Q( m (N)/3. Finally, with the tuning condition (5.47) and 
m = 1/2 we have 

U 2 + 8C(A0 - 45 V / C(iV) (n - ^) + O(6 3 ,5 2 j0N),5aN)X 3/2 (N)) = (5.57) 

with the solution 5 ~ ^J((N) and z — 1. 

Note that one arrives at the same results by considering fcj = 1/4 for i = 1, . . . , d — 1 and 
kd = 1/4 — 1/N, corresponding to momenta p 2 (k) & 2d, for which 4e 2 (/c) ~ 1 [see comment 
after eq. (5.33)]. In fact, in this case, we have 

c(k) = + G(C 3/2 (iV)) (5.58) 

and 

e,(fc) = - [1 - C(N)} . (5.59) 

In particular, for this value of k one can verify that the equation (5.40) has solutions A = 1 — 5 
with 5 finite. Thus, the corresponding relaxation times r do not diverge and CSD is related to 
the solutions of the equation (5.39), yielding z not smaller than 1. 

We checked numerically these results in two, three and four dimensions at (5 — oo for several 
lattice sides, obtaining indeed z — 1 and a computational cost equivalent to that of the standard 
overrelaxation method. [This check was done by updating the components of the f a (k) field in 
the order even-even, odd-even, even-odd and odd-odd.] 

To sum up we can say that with the two matrices of size 4x4 considered in this section, 
we cannot get a dynamic critical exponent z smaller than 1. In particular, we found that the 
characteristic equation of the first 4x4 matrix [see eq. (5.28)] can be factorized in two different 
second-order equations — as predicted in Ref. [10] — giving eigenvalues identical to those found 
when considering a matrix of size 2x2 [see Section 4]. For the second 4x4 matrix [see eq. 
(5.34)] the same factorization leads to two "almost" second-order equations and one can check 
that, compared to the 2x2 case, the extra term (proportional to a/A) does not really modify 
the structure of the equations. Thus, the dynamic critical exponent z must be the same in the 
two cases, i.e. not smaller than 1. 
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5.3 Generalized Stochastic Overrelaxation Algorithm 

One can try to generalize the stochastic overrelaxation algorithm by considering a probability 
distribution different from p m (1 — p). For example, let f(p,m) be the probability of having m 
micro-canonical sweeps [see eq. (4.72)] followed by one standard Los Alamos update. Then one 
has 

oo 

(m) = Yl mf(p,m) = m(p) , (5.60) 

m=0 

where we suppose that Y^m=o f{Pi m ) = 1- (Here p is a tuning parameter or a set of tuning 
parameters.) Note that, if we set f(p,m) = 5 m0 where 5 m o is the Kronecker delta, we obtain the 
Los Alamos algorithm, i.e. we don't do any micro-canonical update and m{p) = 0. 

For this generalized stochastic overrelaxation algorithm we should consider the matrix 

oo 

M(k,p) = f(p,m)M(k,m) (5.61) 

m=0 

instead of the matrix defined in eq. (4.87). It is easy to check that this matrix is given by 

cosJW/2]/ -£i(*(*),p) £ 2 (W,P) \ ( , 62) 

{ ' P) sm[9(k)/2]\-cos[9(k)/2}X 1 (9(k),p) cos \9{k) / 2] E 2 (0(fc),p) J ' l ' ' 



where 



Ei(0(ib),p) = f(p,m) sin [m 9 (k)} (5.63) 



m=0 



^2(9(k),p) = f(P> m ) sin 



m=0 



m + ^ ) f)(k ) 



(5.64) 



The matrix A4(k,p) has eigenvalues and 



mk) ' p) = Zto/2] {cos[^W/2]S 2 (g(fe),p) - U9(k), P )} (5.65) 

= - sm[9(k)/2]E 1 (6(k),p) + cos 2 [^(A;) / 2] S 3 (^(A;),p) , (5.66) 

where we define 

oo 

E 3 {0(k),p) = f^ m ) cos[m9(k)] . (5.67) 

m=0 

For f(p,m) = 5 m o we have J2i(9(k),p) = 0, £ 3 (#(A;),p) = 1 and [using eq. (4.75)] 

\(6(k),p) = cos 2 [9{k) 1 2] = 4c 2 (A;) , (5.68) 

in agreement with the result obtained for the Los Alamos algorithm in Section 4.1. 

In order to study the CSD of this algorithm we should consider the largest eigenvalue of 

the matrix M(k,p). It is obvious that, if 9(k) = [corresponding to c(k) = 1/2], we have 
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£3(0, p) = 1, Si(0,p) = and \(0,p) = 1. At the same time, for a small angle 9{k) = 9 S1 we can 
expand the above expressions in powers of 9 S and obtain 15 

M0*,P) = hi™) + 0(6 3 s ) (5.69) 
M0s, P ) = 1 - |(m 2 ) + 0{6t) (5.70) 

K0 S ,P) = 1 - |((m) + (m 2 ) + i) +O(0t), (5.71) 



where ( to ) and ( m 2 ) are, in general, functions of p. This implies 

r 112 



(to) + 1 1 - X(9 s ,p) 9 2 (m) + (to 2 ) + \' 
From eqs. (4.31) and (4.75) and the relation 9{k) = 9 S we have that 



(5.72) 



1 fi7T 2 

e 2 s *8((N) = - m . (5.73) 

Thus, if ( to ) and ( m 2 ) stay finite we have that r oc iV 2 and z — 2. In order to reduce CSD 
we have to tune p so that both ( m ) and ( m 2 ) go to infinity as powers of 1/9 S . Note that if by 
tuning p we have that ( to ) goes to infinity, then from the inequality ( to ) < ( m 2 ) , which is a 
consequence of the positiveness of the variance a 2 = ( m 2 ) — ( m ) 2 , we get that ( m 2 ) goes to 
infinity too. The same inequality allows, in principle, ( m 2 ) to go to infinity while ( m ) stays 
finite, but this cannot happen if a 2 is finite. Moreover, if a 2 < +00 and ( m ) ~ 9j n we should 
have that ( m 2 ) ~ 6>~ 2n and we obtain 

2 1 2 
r ~ p 1 - < m2 ) ~ • ( 5 - 74 ) 

From eqs. (5.69) and (5.70) and from the fact that T>i(9(k),p) and £3 (#(&), p) are finite, 16 we 
obtain that — at least to order 9 2 — the only tuning we can have is given by 

H - 0^ (m 2 ) ~ 9~ 2 , (5.75) 

which implies [see eq. (4.107)] 

r - 2 ~ ^ dN (5 76) 

r ~ - ~ ^— . (5.76) 

This yields, as expected, z — 1. Thus, with an appropriate choice of the distribution f(p,m) 
and of the tuning, one can only hope to reduce the factor that multiplies l/9 s ~ \fdN/ix in the 
above equation, but there is no way of having z < 1. 
We can verify the tuning relations 

(m 2 ) ~ (m) 2 ~ iV 2 (5.77) 



15 A similar analysis can be done for the case 9(k) = 2tt — 6 S , with 9 S small, corresponding to c(fc) w —1/2. 

16 Actually, from eqs. (5.63) and (5.67), from the relation f(p,m) > and the normalization condition for the 
probability distribution /(p, m) we have that |£i(#(A;),p)| and \T,3(9(k),p)\ are smaller than or equal to 1. 
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for the probability distribution f(p,m) = p m {l — p), which has a tuning condition P ~ 1/N (see 
Section 4.3), and for the probability distribution f(p,m) constant in the interval [1,2m — 1], 
which has a tuning condition m ~ N (see Ref. [9]). In the first case we know that 

(m) = (5.78) 

and one can check that 

<m 2 > = P{1+P) = i— P - (5 79) 

{ ' (1 -pf 2 Pi ■ [ } 

Thus, when P goes to zero as 1/N, relations (5.77) are satisfied. In the second case we have that 
the average value of m is m and one can verify that 

( m 2 ) = - (4m - 1) . (5.80) 
o 

Thus, when m goes to infinity as iV we obtain again relations (5.77). 

Let us notice that the analysis presented in Ref. [10] clearly applies also to the stochastic 
overrelaxation algorithm for any probability distribution f(p, m), as long as the updating matrix 
M(k,p) is 2 x 2. So, the result that z cannot be smaller than 1 for the generalized stochastic 
overrelaxation algorithm given by the matrix (5.62) is not unexpected. However, we believe that 
the previous analysis, and especially the relations in eq. (5.75), clarify how critical slowing-down 
is reduced by this algorithm. 



6 Numerical Results 

In order to check the analytical predictions presented in the previous sections we have done 
numerical tests in two, three and four dimensions. In each case we considered eight different 
lattice sides N, namely N = 16, 32, 48, ... , 128 in two dimensions, N = 8, 16, 24, ... , 64 
in three dimensions and N — 4, 8, 12, . . . , 32 in four dimensions. Also, for all the algorithms 
we have done tests using both the lexicographic and the even/odd update. For the Fourier 
acceleration method only lattice sides that are powers of 2 were considered and we used either 
the whole lattice or even/odd sublattices to implement the Laplacian preconditioning. 

Simulations were done on the PC cluster installed in July 2001 at the IFSC-USP in connection 
with a grant from FAPESP ( "Projeto Jovem Pesquisador" ) . The system has 16 nodes and a server 
with 866 MHz Pentium III CPU and 256/512 MB RAM memory (working at 133 MHz) and is 
operating with Linux Debian. The machines are connected with a 100 Mbps full-duplex network. 
The total computer time used for the tests (including the runs described in Sections 5.2, 6.3, 6.4 
and 7) was equivalent to about 100 days on one node. 

In Ref. [4] we checked the convergence of the gauge-fixing algorithms by considering six 
different quantities. We found that, for each given algorithm, these quantities relax to zero with 
the same speed, i.e. the same relaxation time. Here we consider only two of these six quantities, 
namely 17 

(V^) 2 = E [(V-A c )(x)] 2 , (6.1) 

V x c =l 



17 



In Ref. [4] these two quantities were called, respectively, e-i and e&. 
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which is commonly used in numerical simulations, and 



J <9 



Y d 3 N 



d 



E^E E [QlM-Ql] 2 [Q r 

H=l ° iv c=l x M =l 



(6.2) 



which provides a very sensitive test of the goodness of the gauge fixing [4]. Let us recall that we 
defined the lattice gauge field as [see eq. (2.13)] 



and that [see eq. (2.12)] 



is the lattice divergence of 



(V-A c )(x) = Y: [Afc) - Afa - ej 



A;(x) = -Tr[A,(x)a c ] , 



(6.3) 



(6.4) 



where a c is a Pauli matrix and c = 1,2,3. We also define 



JV 



Q/J, — jy E QfJL^ X ^) 5 



X u =l 



where the quantities 



QIM = E E A » 



// = !,... ,d 



(6.5) 



(6.6) 



(6.7) 



are constant, i.e. independent of x^, if the Landau-gauge-fixing condition is satisfied [4]. Let us 
notice that, at (3 = oo and at a minimum of Eu [g], one has Cf^ = 0. Therefore at (3 = oo the 
quantity Sq should be defined as 



-i d -i 3 JV 



(6.8) 



We used (VA) 2 < 10 -15 as stopping condition for the gauge-fixing algorithms. The quan- 
tity Sq and the minimizing functional Eu [g] have been evaluated only for the final gauge-fixed 
configuration. 

Let us notice that, using eq. (2.18), we can re-write (VA) 2 as 



(va) 2 = y/E E i wC ( x )] 2 ■ 



X c=l 



Also, from (2.6) and (2.8) we can write 

w(x) = J\f(x) w(x) , 



(6.9) 



(6.10) 
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where w(x) is an SU{2) matrix, so that 



1 3 



V X c =l 

At the same time we can re-write the minimizing functional [see eqs. (1.1) and (2.10)] as 



M 9 ] = 1 - ^EE f + ^»V-eJ j = 1 - y w (l ) (6.12) 



and using eq. (6.10) above we have 



^.[^i-^EW, 



3 

t2 



E[^)] 2 - (6.13) 



Then, for (5 = oo and in the limit of large number of gauge-fixing sweeps t, namely by using eq. 
(3.11), we obtain 



(VA) 2 = yEE(w C Wf + ^) (6-14) 

V X c=l 



and 
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M<?1 = ^EE[* c Wf + ^ 2 ). (6-15) 

Thus, with the gauge fixing, (VA) 2 goes to zero and we should find (VA) 2 ^> Sjj [g] in the final 
gauge-fixed configuration. 

For the quantity (VA) 2 (t), in the limit of large number of gauge- fixing sweeps t, we introduce 
a relaxation time r through the relation [see eq. (3.36)] 

(VA) 2 (t) « b exp (- 1 1 t) . (6.16) 

The evaluation of r is done using a chi-squared fit of the function log(VA) 2 (t) . In order to 
get rid of the initial fluctuations, this fit has been done five times using the data corresponding 
to t > n g f(l — l/tfac), where n g f is the total number of sweeps necessary to fix the gauge and 
tf ac = 2,4,8, 16 and 32. In most cases, r increases for increasing tf ac , reaching a plateau. We 
have chosen as the final value for r the second point where the r values become stable within 
errors. 

For the algorithms depending on a parameter, and therefore requiring tuning, we used a 
procedure in three steps in order to find the optimal choice of the parameter, namely the value 
that minimizes the relaxation time r at a fixed lattice side N. This procedure is similar to the 
one described in Ref. [4]. We considered, respectively, 5 configurations in the first step, 10 in the 
second and 20 in the third and final step. 

Our final data for the relaxation time r, the number of gauge-fixing sweeps n g f and the time 
t g f (measured in seconds) necessary to complete the gauge fixing are reported in Tables 1-5 for 
the two-, three- and the four-dimensional cases. When necessary we also report the optimal 
choice for the tuning parameter. For all these quantities we don't show the statistical error since 
it is usually very small, namely less than 1%. Recall that at (3 = oo the link variables U^(x) are 
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set equal to the identity matrix and the initial {g(x)} configuration is chosen randomly. Also, as 
said in the Introduction, in the limit of large number of gauge-fixing sweeps t the configuration 
Ujf\x) = g(x) g~ 1 (x + e IJ ) is driven by the gauge fixing to the vacuum configuration Ujf\x) = 1L, 
losing memory of the initial {g(x)} configuration. This explains why there are usually very small 
fluctuations in the quantities that one is interested in. 

6.1 Critical Exponents and Computational Cost of the Algorithms 

From the data shown in Tables 1-5 one can evaluate the dynamic critical exponents z for the 
five algorithms. In all cases — and in two, three and four dimensions — the results are in 
agreement with our findings in Ref. [4], namely z ~ 2 for the Los Alamos method, z ~ 1 for 
the three improved local algorithms and z xs for the Fourier acceleration method. One can 
also check that the total number of gauge- fixing sweeps n g f grows approximately as N z , with 
the same values of z given above. The time t g f (measured in seconds) necessary to complete the 
gauge fixing grows approximately as N z+d , as expected. The only exception is the Los Alamos 
method, for which we find n g f ~ iV L5 and t g f ~ suggesting that in the initial gauge- 

fixing sweeps this method is more effective than the other local algorithms. This is not surprising, 
since it is well known [2] that the optimal strategy for the overrelaxation algorithm is precisely 
to vary the parameter oo from an initial value 1 (corresponding to the Los Alamos method) to a 
larger asymptotic value u opt . The computational cost of the Fourier acceleration method will be 
discussed in more detail in Section 6.4 below. 

We have also looked at the values of Eg, (VA) 2 and £u[g] in the final gauge-fixed configu- 
rations. We observed that: 

• For all gauge-fixing algorithms and all dimensions considered one finds £q > (VA) 2 > 

£u[g\. 

• These inequalities become stronger as the lattice side N increases. 

• As found in Ref. [4], the Fourier acceleration method is very efficient in relaxing £q and 
in this case one finds £q > (VA) 2 . 

• The ratio Sq/(VA) 2 is usually smaller (or much smaller) for the even/odd update than for 
the lexicographic update. 

• If one uses the quantity £<g to check the convergence of the gauge fixing, then the stochastic 
overrelaxation algorithm is better than the other local algorithms when considering the 
lexicographic update, in agreement with our findings in Ref. [4]. On the contrary, if one 
considers the even/odd update, then the quality of the gauge fixing for the Cornell method 
becomes almost as good as for the stochastic overrelaxation update. 

• For the three improved local algorithms there is in general a gain in computational cost 
when using the even/odd update compared to the lexicographic update. For the Los Alamos 
method and the Fourier acceleration method the situation is reversed. (See next Section 
for a discussion of the Fourier acceleration method.) 

We can conclude by saying that among the local algorithms the best appears to be the Cor- 
nell method with even/odd update. In fact, from the point of view of computational cost, this 
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method is equivalent to the overrelaxation algorithm and almost twice as fast as the stochastic 
overrelaxation algorithm. At the same time, the quality of the gauge fixing, especially when 
considering the relaxation of the quantity Sq, is better than what is obtained with the over- 
relaxation algorithm and almost equivalent to the performance of the stochastic overrelaxation 
algorithm. Let us recall that, in the limit of large number of gauge-fixing sweeps t, the overrelax- 
ation method and the Cornell method coincide [see eqs. (3.17), (3.18) and (3.23)] at leading order 
in e. Therefore, the different performances of the two methods could be related to a difference in 
the behavior for small t. This possibility is discussed in [4, Section 7.1] and we plan to investigate 
it further. 



6.2 Tuning of the Algorithms 

In this section we check numerically the analytic predictions for the tuning of the three improved 
local algorithms and of the Fourier acceleration method obtained in Sections 3 and 4. 

Overrelaxation method: In this case we have the tuning condition [see eqs. (4.51) and 
(4.59)] 

w " = i + L/n (fu7) 



with 



2 7T 

C °P t ~ ^172 • ( 6 - 18 ) 



In order to find the constant C opt one can write 

Copt = N l^Hl , (6.19) 

^opt 

which can be use to fit the numerical data. In this way, considering the optimal choice of 00 
for the three largest lattice sides N, we find for the lexicographic update C opt = 5.0 ± 0.2 in 
two dimensions, 4.26 ± 0.08 in three dimensions and 3.65 ± 0.05 in four dimensions. The same 
fitting procedure gives C opt = 4.01 ± 0.03 in two dimensions, 3.63 ± 0.04 in three dimensions 
and 3.13 ± 0.02 in four dimensions for the even/odd update. Notice that from eq. (6.18) above, 
which is valid for the even/odd update, we have the analytic predictions C op t ~ 4.44 for d = 2, 
Copt ~ 3.63 for d = 3 and C opt ~ 3.14 for d = 4, in good agreement with our numerical results. 

Cornell method: In Ref. [4] we have found the relation 

Uopt = a op t(M) = a opt 2d(l - (£ min )) , (6.20) 

where ( S m i n ) is the average value of the minimizing functional at the minimum. At j3 = 00 
one has ( E m in ) = and the previous relation becomes u opt = 2 d a opt , in agreement with the 
analysis presented in Section 3 [see eqs. (3.17), (3.18) and (3.23)]. From Tables 2 and 3 one can 
check that this relation is very well verified by our data in the two-, three- and four-dimensional 
cases. 

Stochastic overrelaxation method: In Section 4.3 we have seen that the tuning condition 
for the stochastic overrelaxation algorithm is given by [see eqs. (4.93) and (4.106)] 

p = krr (6 - 21) 
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with 



P 



2vr 



(6.22) 



VdN ' 



Moreover, by comparison with the overrelaxation algorithm, one can write p pa uj — 1. It is 
immediate to check that this relation is indeed verified by our data (see Tables 3 and 4). 

Fourier acceleration: We did not discuss the tuning of this algorithm in Ref. [4]. The 
theoretical analysis in Section 3 gives the simple result 



for any dimension d. In Ref. [4] we obtained — in two dimensions, at finite [3 and in the limit 
of large lattice sides N - - the value a opt ~ 1.28, in qualitative agreement with (6.23). The 
agreement is better, as expected, at (3 — oo. In fact, using the lexicographic update we find 
that a opt = 1 for any lattice side iV and dimension d. On the contrary, using the even/odd 
update we have a opt ~ 1.1, with a slow decrease of a opt as N increases. Also note that n g f 
in the lexicographic case is about two times smaller than for the even/odd update. Thus, for 
the Fourier acceleration method, the even/odd update does not help the convergence of the 
algorithm. This result can be understood if one checks the size of the Fourier components of the 
lattice divergence V • A. In particular, one can check that with the even/odd update the slowest 
relaxing mode corresponds to the shortest wavelength. It is this mode that makes the Fourier 
acceleration method perform worse in the even/odd case. This is related to the fact that the 
even/odd update couples the low- and high-frequency modes [15]. 

6.3 Overrelaxation and Stochastic Overrelaxation without Tuning 

In order to check the analytic predictions presented in Sections 4.2 and 4.3, we studied numerically 
the performance of the overrelaxation and of the stochastic overrelaxation algorithms also in the 
case without tuning. To this end, we did tests in the two dimensional case (at (3 — oo with 
lattice sides N = 16, 32, 48, . . . , 128). For both algorithms we considered the two limiting cases 
studied analytically, namely: 

• for the overrelaxation algorithm we used uo = 1.98 (corresponding to the small value Vl pa 
0.01) and uo = 1.3 (corresponding to VL pa 0.54); 

• for the stochastic overrelaxation we set p = 0.96 (namely P pa 0.02) and p = 0.3 (corre- 
sponding to P pa 0.54). 

As a result we got that, for small values of Q (respectively P), r and n g f are indeed constant 
in N. In particular, for the overrelaxation algorithm, we found r pa 50 and n g f of the order 
of 1500 — 1600 using the lexicographic update and r pa 25 and n g f pa 800 using the even/odd 
update. For the stochastic overrelaxation we have r pa 30 and n g f of the order of 1100 using the 
lexicographic update and r pa 25 and n g f pa 900 using the even/odd update. 

For large values of Q (respectively P), the relaxation time r is well fitted (for both algorithms 
and with both types of update) by r pa 0.0136 N z with z pa 2. We have also checked that in this 
case these two algorithms (without tuning) are better than the Los Alamos method, showing a 
relaxation time about two times smaller. 



C^opt 



= 1 



(6.23) 
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6.4 Alternative Implementations of the Fourier Acceleration Method 



In Ref. [17] we introduced a new implementation of the Fourier acceleration method, in which 
the inversion of the Laplacian is done using a multigrid algorithm, avoiding the use of the fast 
Fourier transform. This makes the method more flexible, i.e. it can efficiently work with any 
lattice side N and not only with N equals to a power of 2. Moreover, the new implementation 
is well suited for vector and parallel machines. In particular, we checked that the computational 
cost shows a linear speedup with the number of processors on an APE100 machine for the four- 
dimensional 577(2) case. In that article, we have also implemented a version of the method using 
conjugate gradient instead of multigrid, leading to an algorithm that is efficient at intermediate 
lattice volumes. 

In this Section we want to compare the performance of the Fourier acceleration method based 
on a fast Fourier transform (FFT) to the performance of the alternative implementations based 
on multigrid (MG) or conjugate gradient (CG) algorithms. To this end, let us recall that, even 
though the overhead for the MG or the CG routine is likely to be larger than the one for FFT, one 
can hope to reduce it by exploiting the fact that multigrid and conjugate gradient (as opposed to 
FFT) are iterative methods. In particular, by changing the stopping criterion for the inversion, 
the accuracy of the solution can be suitably varied, while for FFT the accuracy is fixed by the 
precision used in the numerical code. This is important, since the tuning of the parameter a is 
usually done only up to an accuracy of a few percent. Thus, the inversion of the Laplacian most 
likely will not require the high accuracy employed in the FFT case, making possible a substantial 
reduction of the computational cost. We checked in Ref. [17] that this is indeed the case and we 
found that, with an accuracy of about 10 -5 for the inversion, one obtains an algorithm equivalent 
to the original one (based on FFT). For the present paper we checked again this result and found 
a small bug in the code previously used. After correcting it we got that the accuracy necessary 
for the inversion is about 10~ 3 , yielding a substantial gain with respect to the old result and in 
agreement with the intuitive argument above. In order to compare different implementations of 
the Fourier acceleration method we used here the stopping criterion 

— < KT 3 . (6.24) 

where r t is the magnitude of the residual after t iterations. Recall that we want to solve the 
equation 

-Af(i) = (V- A) c (x) , (6.25) 
where c (x) is the desired solution. Then, the residual is defined by 

r c (x) = (V • A) c (x) + A <p c (x) . (6.26) 

In particular, we have tested six different implementations of the Fourier acceleration method 
described below, in addition to the original version (denoted by FFT- FA), based on FFT (working 
in single precision). 

1. MG-FA: The inversion of the Laplacian is done using MG in single precision; as in Ref. [17], 
we used a W-cycle with 2 Gauss-Seidel sweeps before coarsening and 2 after coarsening. 

2. MGCG-FA: Same as the previous algorithm, but with a CG iteration applied on the 
coarsest grid instead of the Gauss-Seidel iteration. This should allow larger coarsest grids, 
which may be useful if one wants to parallelize the code. 
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3. MGCGEO-FA: Same as the previous algorithm, but using a CG iteration on the coarsest 
grid with even/odd preconditioning. 

4. CG-FA: The inversion of the Laplacian is done using CG in single precision and in the 
stopping criterion we compute explicitly the residual defined in eq. (6.26). 

5. CGr-FA: Same as the previous algorithm, but now in the stopping criterion we use the 
magnitude of the residual vector built by the CG method. Note that a CG method usually 
stores three vectors at each step: the approximate solution x, its residual r and a search 
direction p. 

6. CGrEO-FA: Same as the previous algorithm, but now the CG is done with even/odd 
preconditioning. 

We used these six algorithms for numerical tests at (3 = oo in two, three and four dimensions. 
For the tuning parameter a we used the optimal choice obtained for the FFT-FA algorithm (see 
Table 5). We found that the three methods using MG are practically equivalent, with MG-FA 
and MGCG-FA slightly faster than MGCGEO-FA. On the contrary, among the algorithms using 
CG, the last one, namely CGrEO-FA, is always faster than the other two, but still slower than 
MG-FA and MGCG-FA. If we compare these algorithms to the original FFT-FA we see that 
because of the use of FFT the algorithm gets progressively worse as the lattice dimension d 
increases. In particular, we get that in three and in four dimensions the MG-FA algorithm is 
already faster for lattice sizes 32 3 and 16 4 and the gain is larger when considering the even/odd 
update. 18 The CGrEO-FA algorithm is essentially equivalent to FFT-FA in four dimensions for 
lattices 16 4 or larger. 

Finally, we compared the computational cost of the FFT-FA and the MG-FA algorithms with 
the best among the improved local methods, namely the Cornell method. In Fig. 2 we plot the 
CPU time needed to gauge-fix a configuration using these three methods in the four-dimensional 
case using even/odd update. We can see that the MG-FA method is the fastest already for the 
lattice volume 16 4 . This happens in three dimensions at 32 3 . We note that this analysis is very 
machine- and code-dependent, and that the Fourier acceleration methods are particularly well 
suited for the case (3 = oo. As we noted in Ref. [5], the performance of the FA method is poor 
for small (3, reaching z — 1 at f3 — 0. We are currently investigating this matter. 

7 A Gauges 

The analytic study presented in Sections 2, 3 and 4 for the Landau-gauge minimizing functional 
at f3 = oo — namely when all the link variables U^(x) are equal to the identity matrix 1L — can 
be easily extended to the so-called A-gauges, which have been recently used in several analytic 
[18] and numerical articles [11, 19]. 

To this end, let us recall that a general A-gauge can be defined by considering the minimizing 

18 Note that when using the Fourier acceleration method at finite (3 one should avoid updating simultaneously 
all sites of the lattice, since the resulting move in configuration space might be too large and affect the convergence 
of the method. 
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functional 



£u,x [g] = i - 



Tr 



(7.1) 



Clearly, if A M = 1 for all /x we get back the standard Landau-gauge minimizing functional given 
in eq. (1.2). Also, if A« = 1 for i = 1,2, . . .,d — 1, then we can interpolate [18] between the 
Landau and the Coulomb gauge by varying A^ between 1 and 0. 

One can easily redo all the analysis in Sections 2 and 3 and observe that all the formulae are 
still valid if we make the substitutions 



d d 

E - Y.K 

/i=i fi=i 



and 



In particular we have 



d 



EV 



(7.2) 



(7.3) 



h(x, A) = J2 A m [ Un( x ) 9^( x + e /J + u l(x - e M ) g ] (x - e M ) 



and, after setting U^(x) = 1L, 

£ u,\ [g] 



2 V (Ej =1 A,) V 2 ^ 
Also, the Laplacian A becomes a A-Laplacian defined by the relation 

(-A A /)(x) = £ A„ [2/(x) - /(x + e „) - /(* - e„) 

with eigenvalues in momentum space given by 



(7.4) 



(7.5) 



(7.6) 



Pl(k) =4: ]T K sin 2 (7rA; M ) . 

Thus, the Fourier acceleration method is now a A-Laplacian preconditioning, i.e. 

1 



u(x) 



F- 1 



[Pl(k) 



Fw 



(x 



Finally, by using eq. (3.4), we can rewrite the minimizing functional (7.5) as 



£x\f\ 



2V (EUK) 



E /(*)•(- A * 7) (*) +^ 3 )> 



(7.7) 



(7.8) 



(7.9) 
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which implies 



h(x) = 2 £A M 1L - tea- £ A M [ f(x + e„) + f(x - e„) j + 0(e 2 ) 



(7.10) 



w(x) = 2 £A M + K [ 2 /(*) - /(* + e„) - f(x - e„) ] + 0(e 2 ) (7.11) 



v//=l / /i=l 



and 



A/» = 2 (^A.J + 0(e 2 ) 

w(x) = e (-A A /)(x) + 0(e 2 ) . 
Then, the update for the five gauge-fixing algorithms can be written as 
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-t(LosAi.) 

I / (a;) 



with probability 1 — p 



(7.12) 
(7.13) 

(7.14) 

(7.15) 
(7.16) 
(7.17) 

(7.18) 



= (1 - «) /(*) , 

and all the observations reported at the end of Section 3 still apply including the relation 

lo = aAf(x) +C(e 2 ) . (7.19) 

The formulae in Section 4 are also unchanged, with the only exception of eq. (4.29), which 
now becomes 



c(k,X) = 



1 



1 



(X^=i A M ) /i 



£ A M cos ( 2 Ti k h 



1 - 



2 (Ej =1 A,) 



(7.20) 



For the smallest non-zero momentum, in the limit of large lattice side N, this gives 19 



c(k(N),X) 



2Avr 2 



(m=i a,) 



gll - C(#,A)] , 



(7.21) 



19 One can check that this relation holds also for the largest momentum if one considers the absolute value of 
c(k(N),X). 
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where A = min M A M . By using eq. (4.28) we can write 

C(N,X) = (7.22) 

with 

/(A) ee 5^ . (7.23) 
Note that eq. (7.21) is still valid in the case of asymmetric lattices if we set 

& - m „ in % ■ (™) 

It follows that the analysis of CSD and, when necessary, of the tuning of the local algorithms 
considered here is modified in the following way: 

dN 2 

T LosAlamos(X) ~ ~ A /(A) = T^ os Alamos f (A) (7-25) 

"Wr(A) « ^ « //(A) = w //(A) (7.26) 

w(A) « i « ^TV^ = r^y/fiXj. (7.27) 

Clearly, the relation p « cj — 1 is still valid. Also, if we fix the lattice sides iV, the quantity /(A) 
- and therefore the relaxation times of these algorithms — increases (respectively decreases) by 
decreasing (respectively increasing) the value of A. This confirms the results obtained numerically 
in Ref. [11] for the cases Aj = 1 for % — 1, 2, 3, and A 4 = 1 and 0.5. 

In order to verify these results we have done numerical tests in the two-dimensional case 
with Ai = 1 and A 2 = 0.25 for the lattice sides N = 16,32,48, . . . , 128. This choice of A's gives 
H{x) ~ 1.25 and /(A) = 2.5. The data are reported in Tables 6-10. [Again we don't show the 
statistical error since it is usually very small.] By comparison with the data obtained in Landau 
gauge and reported in Tables 1-5, one can easily check the relations for the relaxation time r 
given in eqs. (7.25)-(7.27) above. Also note that eq. (7.19) and the relation p ps uo — 1 are very 
well satisfied by our data. 



8 Conclusions 

We studied numerically and analytically five gauge-fixing algorithms in £77(2) lattice gauge 
theory by considering the case (5 — oo, for Landau gauge and A-gauges. The analysis has been 
done for general dimension d and numerical checks were carried out at d = 2, 3 and 4. Results 
are in agreement with those obtained previously in Landau gauge at finite (5 in two dimensions 
[4]. In fact, we find that the (local) Los Alamos method has dynamic critical exponent z ~ 2, 
the three improved local methods we considered — the overrelaxation method, the stochastic 
overrelaxation method and the so-called Cornell method — have critical exponent z ~ 1, and 
the global method of Fourier acceleration completely eliminates critical slowing-down. 
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As said in the Introduction, if the system does not undergo a phase transition going from 
f3 = to P = oo, then the dynamic critical exponent z should not depend on the constant 
physics, i.e. it should be the same at finite (3 and at f3 — oo. On the contrary, the constant c 
obtained from the fit r = c N z should be different in the two cases and one expects 

c(P = oo) < c( finite (3) . (8.1) 

To make this comparison simpler, we report in Table 11 the values obtained in Ref. [4] of the 
constant c for the five gauge-fixing algorithms at finite (3 and the results of the fits done for the 
same algorithms at (3 — oo. (In both cases we consider the lexicographic update.) From the data 
it is clear that the constant c satisfies very well the above inequality for the five algorithms. 

Our numerical simulations show that the Cornell method with even/odd update is the best 
among the local algorithms. It is very fast and at the same time effective in relaxing the value of 
the quantity Sq. It would be interesting to check if this is true also for finite (3 and for the SU (3) 
case. As already observed in Ref. [4] , among the local algorithms one should choose the stochastic 
overrelaxation method if the lexicographic update is considered. Finally, as expected, the Fourier 
acceleration method is extremely efficient at (3 = oo and we checked that its implementation can 
be improved by inverting the lattice Laplacian using a MG algorithm. 

The theoretical analysis, valid for any dimension d, helped us clarify the tuning of these 
algorithms. In particular, the relations between the parameter uo of the overrelaxation, the 
parameter a of the Cornell method and the parameter p of the stochastic overrelaxation method 
simplify the tuning and confirm nicely the expressions obtained numerically in Ref. [4]. For the 
Fourier acceleration method we found analytically the tuning condition a — 1. This result is 
well verified numerically at (3 = oo and at finite (3 (see Ref. [4]). 

We also studied generalizations of the overrelaxation and of the stochastic overrelaxation 
algorithms. In particular, following a suggestion in [10], we considered explicitly a local algorithm 
(similar to overrelaxation) corresponding to an updating matrix of size larger than 2x2 (i.e. 
4x4). In all cases, we verified that one cannot have a dynamic critical exponent z smaller than 
1 with these local algorithms. 

To sum up, in this work we have done a careful analysis of CSD for the problem of numerical 
gauge fixing in Landau gauge — for the SU(2) group and in d = 2,3 and 4 dimensions — in 
a case that can be studied analytically, i.e. (3 = oo. This study has provided several analytic 
predictions, which we verified numerically. We note that at (3 = oo these results clearly apply also 
to the problem of numerical gauge fixing in Coulomb gauge. We believe that these predictions 
can be very useful in the investigation of more realistic (i.e. finite) values of (3 as well as for f3 — 
and in the extension of this analysis to the general SU (N) case. 
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gure 1: Plot of the eigenvalue A(r, P) [see eq. (4.96)] as a function of r = | | for the 
= 4 and P = 0.2. 
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Figure 2: Plot of the time t g f (in seconds) used to complete the gauge fixing as a function of the 
lattice side N for the Cornell method (+), the FFT-FA algorithm (□) and the MG-FA algorithm 
(O) in four dimensions and considering even/odd update. 
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Table 1: The relaxation time r, the number of sweeps n g f and the time t g f (in seconds) for 
the Los Alamos method for d = 2,3, and 4 using lexicographic (upper rows) or even/odd (lower 
rows) update. 



49 



V — N d 


a 


r 


n 9f 


tgf 


16 2 


0.410 


2.40 


73.5 


0.025 


32 2 


0.445 


4.27 


131.7 


0.175 


48 2 


0.460 


5.85 


186.7 


0.554 


„ , o 

64 2 


0.470 


7.70 


252.9 


1.800 


80 2 


0.475 


9.30 


306.4 


3.690 


96 2 


0.477 


10.48 


342.4 


6.200 


112 2 


0.479 


11.70 


377.6 


9.200 


128 2 


0.482 


13.58 


443.3 


13.800 


n 

8 3 


0.235 


1.21 


38.4 


0.037 


16 3 


0.275 


2.30 


73.2 


0.835 


„ , 'J 

24 3 


0.290 


3.21 


105.0 


4.230 


„ „ o 

32 3 


0.300 


4.26 


141.8 


13.500 


7775 

40 3 


0.305 


5.12 


170.0 


32.500 


48 3 


0.310 


6.33 


210.0 


71.500 


56 3 


0.310 


6.51 


211.0 


115.000 


64 3 


0.312 


7.30 


238.0 


1 AA AAA 

199.000 


4 4 


0.155 


0.56 


20.1 


0.013 


8 4 


0.185 


1.26 


41.0 


0.660 


12 4 


0.200 


1.82 


59.1 


5.120 


16 4 


0.210 


2.38 


79.0 


21.900 


20 4 


0.215 


2.84 


93.0 


62.800 


24 


a on a 

0.220 


o on 

3.39 


1 1 O A 

112.0 


1 rn AAA 

158.000 


28 


0.225 


4.16 


138.0 


348.000 


32* 


0.225 


4.22 


139.0 


622.000 


n ^ „ a 

16 2 


0.405 


0.87 


36.0 


0.013 


32 2 


0.440 


1.89 


62.0 


0.086 


. ') 

48 2 


0.460 


2.88 


97.0 


0.302 


64 2 


0.470 


3.91 


131.0 


0.906 


80 2 


0.474 


4.74 


154.0 


2.250 


96 


A A 7A 

0.479 


5.88 


1 AO A 

193.0 


A A TA 

4.470 


112 2 


0.482 


6.69 


220.0 


6.760 


12© 


A A O A 

0.484 


7.79 


255.4 


1 A 1 AA 

10.100 


n 

8 3 


0.240 


0.55 


22.2 


0.022 


16 3 


0.275 


1.06 


41.0 


0.535 


„ , o 

24 3 


0.290 


1.67 


58.0 


2.860 


32 3 


0.300 


2.12 


78.0 


9.000 


777* 

40 3 


0.305 


3.04 


94.2 


21.400 


48 3 


0.310 


3.23 


114.0 


45.400 


56 3 


0.312 


3.94 


129.1 


83.000 


64 3 


0.315 


4.42 


147.0 


142.000 


4 4 


0.150 


0.29 


14.9 


0.011 


8 4 


0.185 


0.68 


25.0 


0.485 


12 4 


0.200 


0.96 


35.2 


3.660 


16 4 


0.210 


1.30 


46.5 


15.400 


21) 


a on A 

0.220 


1 oo 

1.82 


a c a 
00. U 


C 1 O AA 

51.8UU 


24 4 


0.225 


2.20 


79.0 


131.000 


28 4 


0.225 


2.28 


79.0 


244.000 


32 4 


0.228 


2.55 


89.0 


473.000 



Table 2: The tuning parameter a, the relaxation time r, the number of sweeps n g f and the time 
t g f (in seconds) for the Cornell method for d = 2,3, and 4 using lexicographic (upper rows) or 
even/odd (lower rows) update. 
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Table 3: The tuning parameter u, the relaxation time r, the number of sweeps n g f and the 
time t g f (in seconds) for the overrelaxation method for d = 2,3, and 4 using lexicographic (upper 
rows) or even/odd (lower rows) update. 
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7.41 


271.6 


267.000 


4 4 


0.160 


0.60 


23.2 


0.015 


8 4 


0.400 


1.13 


43.1 


0.824 


12 4 


0.540 


1.66 


63.1 


6.590 


16 4 


0.630 


2.19 


83.8 


27.900 


21) 


0.690 


2.71 


1 AO O 

103.8 


O A C A A 

o4.5UU 


24 4 


0.730 


3.21 


122.0 


206.000 


28 4 


0.760 


3.71 


140.0 


442.000 


32 4 


0.790 


4.26 


162.0 


882.000 



Table 4: The tuning parameter p, the relaxation time r, the number of sweeps n g f and the time 
t g f (in seconds) for the stochastic overrelaxation method for d = 2,3, and 4 using lexicographic 
(upper rows) or even/odd (lower rows) update. 
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V = N d 


a 


T 


n 9f 


tgf 




1.000 


0.06 


5.1 


0.004 


32 2 


1.000 


0.06 


5.4 


0.015 


64 2 


1.000 


0.04 


6.0 


0.092 


128 2 


1.000 


0.05 


6.0 


0.679 


8 3 


1.015 


0.01 


6.9 


0.011 


16 3 


1.000 


0.06 


5.0 


0.118 


32 3 


1.000 


0.06 


5.0 


2.590 


64 3 


1.000 


0.07 


5.1 


30.500 


4 4 


1.000 


0.05 


5.0 


0.005 


8 4 


1.000 


0.05 


5.0 


0.145 


16 4 


1.000 


0.06 


5.0 


6.890 


32 4 


1.000 


0.06 


5.0 


136.000 


16 a 


1.150 


0.27 


11.0 


0.009 


32 2 


1.140 


0.26 


11.0 


0.042 


64 2 


1.130 


0.25 


10.2 


0.232 


128 2 


1.125 


0.25 


10.0 


2.580 


8 3 


1.150 


0.27 


11.0 


0.024 


16 3 


1.140 


0.26 


10.3 


0.324 


32 3 


1.130 


0.25 


10.0 


8.990 


64 a 


1.115 


0.24 


10.0 


109.000 


4 4 


1.155 


0.27 


11.0 


0.013 


8 4 


1.140 


0.26 


10.4 


0.387 


16 4 


1.125 


0.25 


10.0 


24.000 


32 4 


1.105 


0.23 


9.2 


443.000 



Table 5: The tuning parameter a, the relaxation time r, the number of sweeps n g f and the time 
tgf (in seconds) for the Fourier acceleration method for d = 2,3, and 4 using for the Laplacian 
preconditioning the whole lattice (upper rows) or even/odd sub lattices (lower rows). 
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V — N 2 

V — I V 


/ 


n gf 


l af 


16 2 


16.46 


378.1 


0.115 


32 2 


65.13 


1217.7 


1.460 


48 2 


146.19 


2343.3 


7.310 


64 2 


259.68 


3750.3 


23.300 


80 2 


405.58 


5209.9 


61.700 


96 2 


583.91 


7033.1 


119.000 


112 2 


794.66 


8754.2 


198.000 


128 2 


1037.80 


9996.4 


292.000 


16 2 


16.30 


392.3 


0.125 


32 2 


64.93 


1279.7 


1.590 


48 2 


145.99 


2527.7 


7.060 


64 2 


259.47 


4095.9 


26.400 


80 2 


405.37 


5799.6 


80.000 


96 2 


583.69 


7714.4 


162.000 


112 2 


794.44 


9727.4 


285.000 


128 2 


1037.60 


10001.0 


383.000 



Table 6: The relaxation time r, the number of sweeps n g f and the time t g f (in seconds) for the 
Los Alamos method for d = 2 using lexicographic (upper rows) or even/odd (lower rows) update 
with Ai = 1 and A 2 = 0.25. 
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V — N 2 

V — 1 V 


Cc 


q- 




L gf 


16 2 


0.680 


2.71 


87.1 


0.029 


32 2 


0.730 


5.14 


163.6 


0.217 


48 2 


0.750 


7.42 


235.6 


0.708 


64 2 


0.760 


9.61 


298.9 


1.920 


80 2 


0.767 


11.67 


369.9 


4.670 


96 2 


0.772 


13.86 


439.0 


7.940 


112 2 


0.775 


15.57 


484.5 


11.700 


128 2 


0.777 


17.44 


539.4 


16.700 


16 2 


0.685 


1.60 


50.0 


0.018 


32 2 


0.740 


2.95 


101.0 


0.140 


48 2 


0.755 


4.68 


139.4 


0.437 


64 2 


0.767 


5.82 


189.0 


1.280 


80 2 


0.772 


7.22 


225.3 


3.380 


96 2 


0.777 


8.64 


274.0 


6.130 


112 2 


0.780 


10.05 


309.6 


9.550 


128 2 


0.782 


11.41 


352.9 


14.100 



Table 7: The tuning parameter a, the relaxation time r, the number of sweeps n g f and the time 
t g f (in seconds) for the Cornell method for d = 2 using lexicographic (upper rows) or even/odd 
(lower rows) update with Ai = 1 and A 2 = 0.25. 
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V — N 2 

V — 1 V 


( i J 
UJ 


q- 


a af 


L gf 


16 2 


1.700 


2.67 


90.3 


0.033 


32 2 


1.820 


5.08 


165.5 


0.231 


48 2 


1.870 


7.48 


236.7 


0.747 


64 2 


1.900 


9.46 


310.9 


2.310 


80 2 


1.920 


11.77 


391.2 


4.890 


96 2 


1.930 


13.72 


449.2 


8.380 


112 2 


1.940 


15.91 


525.1 


13.400 


128 2 


1.940 


18.69 


531.9 


17.100 


16 2 


1.690 


0.64 


62.6 


0.024 


32 2 


1.840 


2.83 


98.9 


0.145 


48 2 


1.890 


4.35 


145.9 


0.474 


64 2 


1.920 


5.91 


200.3 


1.490 


80 2 


1.940 


7.99 


267.8 


3.980 


96 2 


1.950 


9.74 


321.0 


7.390 


112 2 


1.950 


10.90 


329.1 


10.500 


128 2 


1.960 


12.24 


400.0 


16.400 



Table 8: The tuning parameter u, the relaxation time r, the number of sweeps n g j and the 
time t g f (in seconds) for the overrelaxation method for d — 2 using lexicographic (upper rows) 
or even/odd (lower rows) update with Ai = 1 and A 2 = 0.25. 
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V = N' 2 


V 


T 


n 9f 


tgf 


16 2 


0.680 


3.42 


121.5 


0.044 


32 2 


0.820 


6.99 


239.2 


0.342 


48 2 


0.870 


10.26 


346.1 


1.260 


64 2 


0.900 


13.69 


455.1 


3.230 


80 2 


0.920 


16.44 


569.9 


7.530 


96 2 


0.930 


20.98 


668.9 


12.700 


112 2 


0.940 


23.92 


779.6 


19.800 


128 2 


0.950 


26.08 


923.8 


30.400 


16 2 


0.670 


2.67 


95.2 


0.036 


32 2 


0.820 


5.16 


185.3 


0.278 


48 2 


0.870 


7.79 


277.9 


1.060 


64 2 


0.900 


10.57 


375.7 


3.040 


80 2 


0.920 


12.62 


446.9 


6.910 


96 2 


0.930 


16.16 


525.6 


12.200 


112 2 


0.940 


18.18 


603.0 


19.100 


128 2 


0.950 


19.75 


703.5 


29.300 



Table 9: The tuning parameter p, the relaxation time r, the number of sweeps n g f and the time 
tgf (in seconds) for the stochastic overrelaxation method for d — 2 using lexicographic (upper 
rows) or even/odd (lower rows) update with Ai = 1 and A2 = 0.25. 



V = N' 2 


a 


T 


n 9f 




16 2 


1.000 


0.06 


5.4 


0.004 


32 2 


1.000 


0.05 


6.0 


0.017 


64 2 


1.000 


0.05 


6.0 


0.093 


128 2 


1.000 


0.06 


6.0 


0.748 


16 2 


1.150 


0.27 


11.0 


0.009 


32 2 


1.140 


0.26 


10.2 


0.040 


64 2 


1.130 


0.25 


10.0 


0.230 


128 2 


1.120 


0.24 


10.0 


1.840 



Table 10: The tuning parameter a, the relaxation time r, the number of sweeps n g f and the 
time tgf (in seconds) for the Fourier acceleration method for d = 2 using for the Laplacian 
preconditioning the whole lattice (upper rows) or even/odd sublattices (lower rows) with Ai = 1 
and A 2 = 0.25. 
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Algorithm 


c at finite (3 


c at infinite (3 


Los Alamos 


0.217 ±0.028 


0.026 ±0.000 


Cornell 


0.609 ±0.183 


0.164 ±0.004 


overrelax. 


0.220 ± 0.049 


0.197 ±0.040 


stoc. overr. 


0.300 ±0.047 


0.152 ±0.021 


Fourier 


2.851 ±0.606 


0.042 ± 0.004 



Table 11: The coefficients c obtained by fits of the relaxation time r for the five gauge-fixing 
algorithms in two dimensions (with lexicographic update) at finite (5 (from Ref. [4]) and at infinite 
P. 



58 



